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ABSTRACT 

We present a new multi-dimensional radiation-hydrodynamics code for massive stel¬ 
lar core-collapse in full general relativity (GR). Employing an Ml analytical closure 
scheme, we solve spectral neutrino transport of the radiation energy and momen¬ 
tum based on a truncated moment formalism. Regarding neutrino opacities, we take 
into account a baseline set in state-of-the-art simulations, in which inelastic neutrino- 
electron scattering, thermal neutrino production via pair annihilation and nucleon- 
nucleon bremsstrahlung are included. While the Einstein field equations and the spatial 
advection terms in the radiation-hydrodynamics equations are evolved explicitly, the 
source terms due to neutrino-matter interactions and energy shift in the radiation mo¬ 
ment equations are integrated implicitly by an iteration method. To verify our code, we 
first perform a series of standard radiation tests with analytical solutions that include 
the check of gravitational redshift and Doppler shift. A good agreement in these tests 
supports the reliability of the GR multi-energy neutrino transport scheme. We then 
conduct several test simulations of core-collapse, bounce, and shock-stall of a 15 Mq 
star in the Cartesian coordinates and make a detailed comparison with published re¬ 
sults. Our code performs quite well to reproduce the results of full-Boltzmann neutrino 
transport especially before bounce. In the postbounce phase, our code basically per¬ 
forms well, however, there are several differences that are most likely to come from the 
insufficient spatial resolution in our current 3D-GR models. For clarifying the resolu¬ 
tion dependence and extending the code comparison in the late postbounce phase, we 
discuss that next-generation Exaflops-class supercomputers are at least needed. 

Subject headings: hydrodynamics — methods: numerical — neutrinos — radiation: 
dynamics — supernovae: general 
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1. Introduction 

Neutrino transport is an essential ingredient for numerical simulations to clarify the the- 


e.g., 

Foglizzo et al. 

(2015 

1; 

(2012b 

) for recent reviews) 


Foglizzo et al.1 (|2015l ) ; iMezzacappa et al.1 (120141 1 ; iBurrowa (120131 1 ; iJankal (120121 1 ; iKotake et al. 
H)for recent reviews). In the collapsing iron core, neutrinos play crucial roles in every phase; 


deleptonization in the core d etermines the time of bounce and the mass of the proto-neutron star 
(PNS) (e.g., Langanke et al. (120031 )): the gigantic internal energy tapped in the PNS is almost 
completely carried away by neutrinos, a tiny fraction of which contributes to the heating of the 
postshock material, leading to the onset of core-collapse supernova e (CCSNe) in the context of the 
neutrino-heating mechanism (jWilsonl Il985l : iBethe &; Wilsonl Il985l l . Since these SN neutrinos are 
generally not in both thermal and chemical equilibrium, the distribution of neutrinos in the phase 
space should be computationally determined. This is a seven dimensional (7D) problem; (3D+1D 
in space-time and 3D in momentum space), the reason why CCSN simulations are considered as 
one of the most challenging subjects in computational astrophysics. 


ri 
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. n>etne 1 

17)) and 

iyyt j 

the st 

nerant et ai. (iyy4 

anding-accretion-sh< 

1 : Burrows et ai. 

3ck instabil- 

ity (SASI, e.g., 

Blondin et al. 

2003); 

Foglizzo et al. ( 

2006,2007) 

; Ohnishi et al. (2006): 

Blondin Sz Mezzacappa 

(2007) 

Iwakami et al. (20081. 2 

009);] 

Fernandez & Thompson ( 

2009) 

Guilet & Foghzzo ( 

2012): Hanke et al. 


( 2012 ); 


Foghzzo et al 


( 20121 ): Couch (20131): [Fernandez et al. 


( 20141 11. the postbounce cores are es¬ 


sentially of multi-dimensional (multi-D) natureQ Due to the high compactness of the PNS, these 
multi-D fluid motions are all under the influence of the general relativistic (GR) grav ity, the consid¬ 


eratio n of which u s ed to be standard in the pioneering era of CCSN sim ulations fe.g..lMav fc White 
(1966); Schwartz ( 1967 )1. In rapidly rotating supernova cores (e.g., Wooslev fc Bloom ( 200(tI )1. 


spherical (e.g.. 

Ardelian et al. 

(2000); 

Kotake et al. 

2004, 

2006b); 

Burrows et al. 

(2007): iMasada et al. 

(2Q12I, 

2014): 

Sawai et al. 

Iwakami et al. 

(2014)). All in all, in order to have the final word 


Obergaulingeret al. ( 

2006, 

2014) 

(2013); 

Nakamura et al. 

(2014al) 


titatively, one needs to deal with the 7D Boltzmann neutrino transport simulations in full GR 
MHD. Unfortunately, however, this still remains to be computationally very demanding even us¬ 


ing exa -scale computing platforms in the next decade(s) to come (see discussions in iKotake et al 
(l2012al ffl 


Since the late 1990s, the ultimate spherically-symmetric (ID) simulations wh ere the GR B oltz¬ 


mann trans port is coupled to 1 D-GR hydr odynamics have been made feasi bl e by 


(I2OOII I and ISumivoshi et al.l (120051 1 (see IMezzacappa fe Matzneri ()1989l l: IMezzacappa &; Bruenn 


Liebendorfer et al 


1 Recently multi-dimensionalities in the precollap s e cor e (e.g., iMeakin et al.1 1 2011 11 are also attracting much 
attention (e.g., Couch fc Ott ( 2f)l.lh : Miiller fc .Tanka (2014)). 

2 We here mean the feasibility of 7D-GR Boltzmann neutrino transport simulation following ~1 s after bounce 
with sufficient numerical resolutions in both space and momentum space. 
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( 1993b 30); Yamada ( 1997 ): Yamada et al. (1999); Bruenn et al. ( 2001 '): Liebendorfer et al. ( 2001 


2004 ) for the code developments). Unfortunately, however, these full-fledged ID sim ul ation s fail 
to produce explosions except for super-AGB st ars at the low-mas s end (jKitaura et alJ 120061 ). In 


the context of the full Boltzmann calculations, Liyne et al. (2004J) were the first to perform two- 


dimensional (2D) multi-angle (i.e., assuming axisymmetry in both space and momentum space) 
neutrino h ydrodynamics sim ulations using the discrete ordinate (S n ) method. Then it was demon¬ 


strated by lOtt et al.l (2008) that the multi-angle transport is really important especially when the 


neutrino radiation fie ld deviates significa ntly from spherical symmetry such as in the rapidly ro¬ 


tating cores (see also |BrandUet_ah| (201 If)). Going beyond the assumption of axisymmetry in the 


multi-angle transport, ISumivoshi Sz Yamadal (120121 1 were the first to develop the fully mul ti-angle 


Boltzma nn transport code and then apply it for static backgrounds (ISumivoshi et al.ll2014l l. More 


recently, [Nagakura e t al. (2014) extended the code to include special relativistic (SR) corrections 


and showed the ability of the code by performing ID core-collapse simulation of a 15 Mq model. Al- 


(20141), and lPere s et al. (20 141) 


Cardall et al. 

(2013b 

a); 

Shibata et al. 


At present, direct discretization of the Boltzmann transport equation fully into the neutrino 
angle and energy (such as by the S n method mentioned above) is still computationally very expen¬ 
sive. An approximation often made in previous works is to remove the full angular dependence of 
the Boltzmann equation by expanding the neutrino distribution function as a series of moments. 
The simplest version, in which one closes the moment expansi on after the zeroth ang ular moment, 


is multi-group flux l i mited diffusion (MGFLD) scheme (e. g .. iBruennl (119850: 


Livne et 


al 


(.20041): 

Kotake et al.l (|2006al ): ISwestv fc Mvral (|2009h : IZhang et alJ (|201Bh : iBruenn et al.l (20131 )). In FLD 


schemes, the basic equation is a diffusion equation for the neutrino e nergy densi t y. In so lving it, an 
appro p riate flux-lim iter should be employed (e.g., iMinerbol (|1978l ): IPomraningl ()198ll ): iLevermore 
(| 19841 ): I.Tankal (119921 11 to ensure the causality of the radiation field in the transparent regions where 
the diffusion approxi mati on b reaks down. The isotropic diffusion source approximation (IDSA) 


scheme developed by IL iebend drfer et ah (2009) is basically positioned at a similar approximation 


level compared to the MGFLD scheme. In the IDSA, the neutrino distribution function is divided 
into two components, the one which is trapped with matter and has isotropic distribution function 
and the other in the free-streaming limit, each of which is solved independently, while satisfying 


been extensively employed in both 2D 

(Suwa et al. 

2010 

2011, 

2013, 

2014: 

Nakamura et al. 

2014b) 

and 3D simulations ( 

Takiwaki et al. 

2012|, 

2014 

). One can also truncate the angular moment at 


the second order and transport the zeroth and first order angular moment. In this case, higher or 
equal to the second order moment needs to be determ ined independentl y to close the s e t of t he 


transport equations. In the Ml moment scheme (e.g., iPons et al.l (j2000l l: IShibata et al.l (201ll ll 


one applies an ana lytic formula for the closu re relation (see examples app lied in post-Newtonian 


MHD simulatio ns (lObergaulinger et al.l 1201411 and GR simulations in ID (lO’Connor &; Ottl 120131 : 


Q’Connorll2014l ) and in 3D (iKuroda et al.ll2012l . 120141 11. In contrast, one can self-consistently de 
































































































































































- 4 - 


termine the closure relation by the variable Eddington factor (VEF) method (e.g., iMtiller et al 
(120101 1 and see references therein). In these cases, a model Boltzmann equation is integrated to 
iteratively obtain the solution up to the higher moments (i.e., the Eddington tensor) until the 
system is converged. Currently, the state of the art in multi-D simulations of CCSNe is defined by 
multi-group (sp ectral) neutrino hydrodynamics simu lations. More severe approximations include 


gray transport (Fryer e 

Janka Sz Muller i 199( k 


al 


1999; 


Scheck et al 


Ruffert et al. (1991): 


200611 or the light-bulb and l eakage schemes (e. 


Murphy fc Burrowsl ( 2008 ): lO’Counor fc Ottl ( 201ll ): IPerego et al.1 ( 20141 )). which have been often 


Rosswog fc Liebendorferl (200.8 ): iKotake et al.l (200 


employed in many recent studies of multi-D instabilities and the MHD mechanism of CCSNe. 

In addition to the multi-D and multi-angle/truncated neutrino transport, the accurate treat¬ 
ment of GR is highly ranked among the to-do lists towards the ultimate CCSN simulations. In 
most of previous multi-D models with multi-group neutrino transport, GR effects are attempted 
to be modelled by using a modified gravitational potential that takes into account a ID GR ef¬ 


fect bv replacing the monopo l e term of Newtonian gravity with the TOY potential (IBuras et al 


2006b c 

: Marek Y Janka 

2009; 

Bruenn et al. 

200d; 

Hanke et al 

20131) 

GR core-collapse simulations in 2D (e.g., 

Dimmelmeier et al. 

20021) 


3) • WhiD t here are a nu mber o f 


(2012b)) and in 3D (e.g., 

Shibata & Sekiguchi ( 

2005 

); 

Ott et al. ( 

2012) 

Kuroda et al. 


(20121,120141); IMosta et a l.l (2014)), many of them especially in 3D have been made, for the sake of 


(Liebendorfer et al. 

2005 

) 

study 

Kuroda et al 

2012) 


20051) or by the neutrino leakage scheme (jSekiguch i I2010I ). In our previous 


with the gray Ml scheme. We demonstrated that due to deeper gravitational well of GR the neu¬ 
trino luminosity and the average neutrino energy in the postbounce phase increase when switching 
from SR to GR hydrodynamics. Since the neutrino heating rates in the postshock regions are 
sensitively affected by the emergent neutrino spectra, the GR effect whether it will or will not help 
the onset of neutrino-driven explosions needs to be investigated by multi-D GR simulations with 
more sophisticated neutrino transport scheme^ 

In this paper, we present a new 3D-GR radiation hydrodynamics code that is meant to apply 
for stellar core-collapse simulations. The spacetime treatment is based on the Arnowitt-Deser- 
Misn er 3+1 formalism and we employ the Baumgarte-Shapir o-Shibata-Nakamura (BSSN) formal¬ 
ism ([Shibata fc Nakamuralll995l : iBaumgarte fc Shapiro! 119991 ) to evolve the metric variables. The 
full GR radiation-hydrodynamics equations are evolved in a conservative form, in which we solve 
the energy-dependent set of radiation moments up to the first order with the Ml moment scheme. 
Tliis part is based on the partial i mplementation of the Thorne’s moment formalism ([Thorne 
19811 ). which is extended by Shibata et al. (201l|) in a more suitable manner applicable to the neu- 


3 It is worth mentioning that 2D GR models with the VEF method t end to explode more easily than the corre¬ 
sponding 2D Newtonian models with and without the GR correction (e.g. Muller et ah 201CI . 2012bl : 


Muller &: Janka 


2014j). This may support the speculation that GR is helpful for the working of the neutrino mechanism in multi-D 


simulations. 
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trino transport problem. Regarding the neut r ino-m atter interaction terms , we e mploy a baseline 
set of weak interactions as given in iBruennl (|1985l ) and iRampp fc Jankal (120021 ) where nucleon- 
nucleon bremsstrahlung is additionally taken into account. Our newly developed code is designed 
to evolve the Einstein field equation together with the GR radiation hydrodynamic equations in 
a self-consistent manner while satisfying the Hamiltonian and momentum constraints. A nested 
structure embedded in the 3D Cartesian computational domain enables us to f ollow the dynamics 
starting from the onset of gravitational collapse of a 15 Mq star (IWooslev fc Weaverl 1995|), through 
bounce, up to about ~50 ms postbounce. Since, it is still computationally too expensive to follow 
long-term evolution in f ull 3D until the neu t rino-driven explosion takes place (e.g., at the earliest 
~ 200 ms after bounce ( Brnerm et a,l. 2009 : Marek &; Janka 2009 ) or ~ 500 ms in 2D GR calcu¬ 
lation (M uller &; Jankal 201^)), we mainly focus on detailed comparisons between our pseudo-ID 
neutrino profiles computed in the 3D Cartesian coordinates and previous ID results to check the 
validity of our new code in the early postbounce phase. 


This paper is organized as follows: In section [2j after we shortly introduce the formulation 
of the GR transport scheme, we describe the governing equations of hydrodynamics and neutrino 
transport in detail. Some practical implementation schemes how to satisfy important conservative 
quantities such as lepton number, energy, and momentum are given in section [3j The main results 
and detailed comparisons with previous studies are presented in section [5l Note that geometrized 
unit system is used in sections [ 2 ] and El i.e. the speed of light, the gravitational constant and the 
Planck constant are set to unity: c = G = h = 1, and cgs unit is used in section [5j Greek indices 
run from 0 to 3 and Latin indices do from 1 to 3. 


2. Formalism 


This section starts with a brief summary of the basic equations and the numerical schemes of 
GR radiation hydrodynamics. 


Following our previous work (jKuroda e t al . 120121 ). our code consists of the following three 
parts, where the evolution equations of metric, hydrodynamics, and neutrino radiation are solved, 
respectively. Each of them is solved in an operator-splitting manner, but the system evolves self- 
consistently as a whole satisfying the Hamiltonian and momentum constraints. Regarding the 
metric evolution, the spatial metric 'jij (in the standard (3+1) form: ds 2 = — o?dt 2 + 7 ij(dx l + 
/ 3 l dt)(dx :) + j3 J dt), with a and f3 l being the lapse and shift, respectively) and its extrinsic curvature 
K, , a re e volved using the BSSN for mulation (jShibata &; Nakamural Il995l : iBaumgarte fc Shapiro 


19991 ) (see Kuroda et al.1 (120121 . 120141 ) for more details). 
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2.1. Radiation Hydrodynamics 


Major differences compared to our previous code (IKuroda et al. 120121 ) are, on the one hand 
we evolved energy integrated ( u gray’ 1 ) neutrino radiation held with an approximate description 
of neutrino-matter interaction based on the neutrino leakage scheme, on the other hand we now 
solve the spectral neutrino transport where the s ource terms ar e treated self-consistently following 
a standard procedure of the Ml closure scheme (jShibata et al. 2011). For convenie nce, we b riefly 


summ arize the basic equations of our newly developed code in the following (see IShibata et al 
(201 lj) and Car dall et al. ('12013 a) for the complete derivation). 


The total stress-energy tensor is expressed as 


nQ;/3 
' (total) 


nOLfH 
' (fluid) 


J de E 


r a/3 

-(*v0’ 


( 1 ) 




where and stress-energy tensor of fluid and energy-dependent neutrino radiation 

held, respectively. Note in the above equation, summation is taken for all species of neutrinos 
(v e ,V ei v x ) with v x representing heavy-lepton neutrinos (i.e. v jt . v T and their anti-particles), and 
e represents neutrino energy measured in the comoving frame with the fluid. For simplicity, the 
neutrino havor index v is omitted below. 


With introducing radiation energy (E^), radiation hux (F^) and radiation pressure (P^p, 
measured by an Eulerian observer or (J( e ), and Tp) measured in a comoving frame, can 
be written in covariant form as 


= E (s)" V + Ffan" + + Pfi, 


flU 


= J {e) v?u v + Hfau" + Ufa? + L {e) 


(2) 

(3) 


In the above equations, n^ = (1/a, — (3 k /a) is a unit vector orthogonal t o the spacel ike h ypersu r face 
and is the four velocity of fluid. In the truncated moment formalism (Thorna ll981 : Shibata ct al 


2011 ), one evolves radiation energy (E^) and radiation hux (P^) in a conservative form and 
is determined by an analytic closure relation (e.g., Eq.([ 6 ])). The evolution equations for and 
FT are given by 


&ty/l E {e) + diy/liaFfo - f3 l E {e) ) + y/^ad^eM^n^) = 
Vl( aP le) K ij ~ F (e) 9 i a ~ aS (e) n m)> 


(4) 


and 

dty/lF (e) . + djypyiaP^ - P J F (e) .) - ^ad £ (sMf E) 7 ^) = 

y/l l-E( e )dia + F( e)j di/3 J + (a/2)P^ ) d i -/ jk + aSfa 7 ^], (5) 

respectively. Here 7 is the determinant of the three metric 7 = det( 7 jj) and is the source term 
for neutrino matter interactions (see appendix [A] for the currently implemented processes). Aph is 
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defined b y gUn where denotes the third rank moment of neutrino distribution 

function (jShibata et al. 0201 1) for the explicit expression). 


By adopting the Ml closure scheme, the radiation pressure can be expressed as 

3x(e) - 1 


_ 

Me) “ 


. pij I ^ Mp) pij 

thin(e) ' 2 thick(e) > 


( 6 ) 


where X( £ ) represents the variable Eddington factor, and P^ lick p) corresponds to the radi- 


and Hi' , ( 

Shibata et a . 

2011 ) 

- M7 - 

Obersaulineer h Janka 

( 201 l|) 


201 1). Following Minerbo (1978); Cernohorskv fc Bludma.nl (1 994 ) and 


X(e) 


5 + SF ( % - 2Ff c) + 6 F* t) 

15 


where 


%) - 


/2 

Me) 


(7) 


( 8 ) 


In Eq.©, = g lw + u^u v is the projection operator. As we will discuss later, by the definition of 

P( e ) in Eq.(|5]). one can appropriately reproduce several important neutrino behaviours, for example, 
neutrino trapping in the rapidly collapsing opaque core. We iteratively solve the simultaneous 
equations J7HH1) to find the converged solution of X(e)- 


The hydrodynamic equations are written in a conservative form as, 


d t p* + di(p*v l ) = 0 , 
dty/xSi + djyfy (S iV i + aPb {) = 


(9) 


—V7 


S 0 dia - S k di/3 k - 2 'aS%di<j> + ae A<t> {S jk - P'y jk )da jk /2 + a J dteS ^ 7 ^ 


dty/jr + d iy /x(rv l + P{v l + (3 1 )) = 


V7 


aKSlt 3 + ae-^iSij - Pj^Av - S l D i a + a j deS {* 


£ ) n M 


dt(p*Y e ) + dfaYeV*) = Viam u J *(5^ - 


( 10 ) 

( 11 ) 

( 12 ) 


where p* = p^fSjW, Si = phWut, Sij = phutUj + Pjij, S k = 7 U , So = phW 2 — P and 
<p = log( 7 )/ 12 . p is the rest mass density, W is the Lorentz factor, h = 1 + e + P/p is the specific 
enthalpy, v l = u l /u t , t = So — pW, Y e = n e /rib is the electron fraction (nx is the number density 
of X), e and P are the specific internal energy and pressure of matter, respectively and m u is 
the atomic mass unit. P(p,s,Y e ) and e(p,s,Y e ) are given by an equation of state (EOS) with s 
denoting the entropy per baryon. We employ an EOS by Lattimer & Swesty (1991) (LS220, see 
section section l5Tl for more detail). In the right hand side of Eg. dill) . D l represents the covariant 
derivative with respect to the three metric 7 ^. 






























2.2. Conservation of Energy and Lepton Number 


As explained in Section f2.11 the formalism of our code that treats the radiation-hydrodynamics 
equations in a conservative form is suitable to satisfy the energy conservation of the total system 
(neutrinos and matters, see also lKuroda &: Umedal (|2010l l: iKuroda et al.1 (|2012l ) for more details). 
Let us first show how the energy conservation law is obtained in our code. 


To focus only on the energy exchange between the matter and neutrino radiation field, we 
omit the gravitational source term in the following discussion. Then, the equations of energy 
conservation of matter and neutrinos (e.g., Eqns. (J4]) and (1111) 1 become 


dt\fl T + di^(TV z +P(v l + /T)) = J deViaS^n^, (13) 

<W l E {e) + diVl( aF {e) - P E (e)) +/ya3 E (^ ) n /1 ) = ~Vl aS ( £ ) n n- ( 14 ) 


From the above two equations, one can readily see that the total energy (sum of matter and 
neutrinos) contained in the computational domain, E um = f dx 3 ^/j(r + f deE ( £ )), is conserved in 
our basic equations as long as there is no net energy flux through the numerical and momentum 
space boundaries (i.e. f de de^eM^n^) = 0). 


The lepton number conservation needs to be satisfied with good accuracy because it determines 
the PNS mass and the postbounce supernova dynamics. We here explain how we treat it in our 
code. As for the electron and neutrino number conservation, the basic equations are given by 


d t 


P*Y e 


TO,, 


di 


p*Y e v l 




de 


= Via / —{S^ £) - S^ e) )u„, (15) 


dt{ql e) Via) + di{q\ e) Via) - Viad £ [eq ^J V 0 u a = 

Q(e) = ~ £ ~ lT (e) U 0 = ( J (e) ua 

n a P — p~ 1 T 1 ^h a — (T-/ a iV - 4 - P a P\ 
q {e) ~ £ 1 (e) n i ~ y H (e) U + L '(e)>' 


Via 


where 


(16) 

(17) 

(18) 


with VJ,'H a , C a P) = e _1 ( J, H a , L a P). T he conservat i on eq uation for neutrino number (1161) cor¬ 
responds to Eq.(3.22) (divided by s ) in [Shibata et ah (2 01 1). Since the neutrino number density 
measured by an Eulerian observer is expressed as 


Euler = J deq°( Uj£) Via, 

the equation of total lepton number conservation becomes 

dt{p*Yi) + di (p*Y e v z + m u Via J de q\ u ^ e) - q\^ e) 


= 0 . 


(19) 


( 20 ) 
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Here the total lepton fraction Yi is defined by 


Yi 


Ye + Y Vr - YL 


Y e + 


Ve 

rn u 
pu l 



v e ,e) ^(V e ,£))' 


( 21 ) 


From Eg. (1201) . one can readily see that the total lepton number is conserved irrespective of the 
included neutrino matter interaction processes in case that there is no net flux through the numerical 
and energy space boundaries. The distribution of Yi into the each component (e.g., Y e ,Y Ue ,Yp e ) is 
determined by the details of the implemented microphysics, which should be checked carefully and 
will be reported in Sec. [5j 


2.2.1. Neutrino Number Transport in the Diffusion Limit 


In the collapsing iron core, it is well known that the central core becomes opaque to neutrinos 
drug mainly to scattering off heavy nuclei when the central density exceeds rsj 10 11 12 g cm 3 (jSatol 
19751 1 and neutrinos are trapped with matter afterward. In the diffusion limit at large neutrino 
opacities, the trapped neutrinos move with the matter velocity Vi for an Eulerian observer. Thus 
their advection equation can be described with the same form of Y e as 

dtp*Y u + dipfY v v l = (Source terms). (22) 


Because the source terms in the above equation amount equal to the negative value of the source 
terms in the electron number conservation equation, the core lepton number is conserved in a good 
accuracy until it gradually decreases by diffusion in the PNS cooling phase. Since the central core 
mass depends on the core lepton number, the CCSN simulation should be able to capture this 
important phenomena appropriately. 

In our formalism, however, this is not a trivial problem because we solve the energy-momentum 
conservation equations Q-dS]) and not the neutrino number conservation equation (|16l) . In this 
section we check whether our basic equations can reproduce the neutrino diffusion limit adequately, 
i.e., they reach asymptotically to Eo. ([22j) . 

For simplicity, we assume in the following that the spacetime is flat and the typical velocity 
of the matter field (v) is much smaller than the speed of light (slow motion limit; neglecting terms 
higher than the second order with respect to (v/c)). 

Let us first check whether Eg. (1161) can satisfy the local neutrino number conservation in the 
trapped region. In this limit, neutrino number density at each energy bin q° and its flux q l approach 

n rr t 1 - 1/0 rr . 

q = JU + H ~ J + —yy~ ~ J , 

q i = Ju { + U l - JN + W ~ gV + W/J). 


(23) 

(24) 
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From these relations, it is obvious that the neutrino number density at each energy bin q° is 
transferred with the matter velocity v l plus the diffusion velocity T~L l / J and the equation of the 
total lepton number (Eq. (1201) ) in the slow motion limit becomes 


dt(p*Yi ) + di ( p*Y e v l + m u y/^a / de ^ 




v l ) = d t (p*Yi) +di(p*Yiv z ) 
= 0 , 


(25) 


demonstrating that Eg. (1201) satisfies the local lepton number conservation in the trapped region. 

Next, we take the diffusion limit of Eq.Q. In this limit, H l /J should approach 0 (i.e., the 
radiation flux (DP) in the comoving frame vanishes). From this, the following relation can be 
derived, 


H l - Ev i + F i — P ij Vj ~ 0 


F 1 ~ -Ev\ 
3 


(26) 


where we take a simple closure relation P lJ = -jin the diffusion limit. Inserting this into the 
left hand side of Eq. (J4|) , one can get 

d t E + d t F l + d e (-eP^d iVj ) ~ d t E + d t - | &d iVj - e (d £ P lj ) d iVj (27) 

~ dtE + di - di - e (d £ pd) d iVj (28) 

- 8 t E + di{Ev i ) - e ( d £ p ij ) d,v r (29) 


Moving from the right hand side of Eq. (1271) to that of Eq. (|28|) , we assumed that E has almost no 
spatial gradient well below the neutrino sphere s (at hig h opac ities) in the prebounce core. This is 
satisfied quite well as shown in the literature (jBruennl (j 19851)), in which a nearly flat Y Ve profile 
is shown in their standard models within the central core with mass ~0.5-0.6 Mq after the central 
density exceeds ~ 5 x 10 13 g cm~ 3 . Note that, in Eq. (1281) . the third term — di(Ev l /3) which is 
originally a part of the advection terms in the energy space d £ (—eP t:, diVj) balances with a part of 
the spatial advection term <9j(4£V/3) and they lead to the second term in Eq. ()29l) . From this, one 
can clearly see how the apparent advection speed of E at each energy bin e approaches the matter 
velocity v l in the diffusion limit. Then, assuming that the neutrino number density at each energy 
bin can be approximately expressed as 


n 

q rsj 


J ~ £ - 2Fvi ~ £, 


(30) 


in the slow motion limit, when we divide Eq. (I29j) by e and integrate it in energy space, it can be 
summarized as 

J de [d t q° + di(q°v l ) - ( d £ P lj ) diVj\ = d t n u + din u v\ (31) 

From this, one can clearly see that the neutrino number density n u (= pY u ) is transported with the 
same matter velocity v l in the diffusion limit. 
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3. Numerical Method 


In this section we describe how to evolve the radiation-hydrodynamics variable^. As we 
explained in the previous section, we solve Eqs. (HJ), (JSJ) and (f9])- (fT2l) as our basic equations which 
are collectively expressed as 


dtQ, + S a( j v s + S a dv,e + S grv + S ism — 0, 


(32) 


where Q denotes conservative variables 


Q 


P* 

v / 7 r 

p*Y e 

- 


(33) 


In Ea. ([32|) . S a d v , S ) S a d v ,e) Sgrv and denote advection term in space, advection term in momen¬ 
tum space, gravitational source and neutrino-matter interaction term, respectively. Throughout 
this paper other than Appendix [0 we divide this equation into the following two parts which are 
expressed in the finite difference expression as, 


for the explicit part and 


Q* - Q n 

At 


I c n 
i °adv,s 


+ s 


grv 


= 0 , 


(34) 


Q«+l _ Q* 

At 


+ S adv , e n+1 + S 


n +1 _ 


= 0 , 


(35) 


for the implicit part, respectively. In Eas. (l34l) - (l35l) . At is the time step size between n-th and 
n + 1-th time steps and the upper indices “n” represents n-th time step. Variables with denote 
the time updated variables during an operator-splitting procedure. 


In Eq. ((34j) . S a d v , s n and S grv n represent the terms with respect to advection in space and 
gravitational fields at n-th time step, both of which are added first in an explicit manner to obtain 
conservative variables at a middle time step Q*. Next in Eq. (1351) . the rest of terms, advection 
in energy space (S a d v , e ) and neutrino-matter interaction terms (S^m) at (n + l)-th time step are 
added to Q* in order to find the converged solution of Q n+1 by an iterative method. The reason 
why we separate source terms into the two parts, explicit and implicit ones, is that the typical 
time step size, At ~ 4 x 10 _7 s, which is determined by the speed of light c, typical minimum grid 
width in our calculation Ax ~ 500m and the Courant-Friedrichs-Lewy number CFL, e.g., 0.25, is 


4 We omit our numerical method to evolve the space-time variables that is essentially the same as in lKuroda et al 


( 2012h . 
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sufRciently short for the advection term in space S a( j v , s an d the gravitational source term S grv as 
well as for all the geometrical variables by an explicit update. However it is too long to follow, 
e.g., the weak-interaction term, which has significantly shorter time scale (< 10 _9 s). We thus need 
to treat these terms in an implicit way through Eg . (135 1) to ensure a numerical convergence and 
stability. 

Here, we should comment on the treatment of the advection terms in energy space S a d v ,e- As 
we mentioned above, we solve them implicitly in time as Eg. (1351) in our main results. It means that 
there is a time gap At between moment of evaluation for advection terms in real space S a d v ,s n and 
that in energy space S a d V e" +1 - Although we can generally obtain consistent results with previous 
studies regardless of the time gap, it is noted that the treatment of the energy advection either by 
the implicit or explicit scheme leads to visible changes in the postbounce features. We will discuss 
this point further in Appendix [Bj 

In the following sections, we are going to describe how to evaluate the advection terms in space 
(Sec. 13.11) and in energy space (Sec. I3.2[) . and then move on to describe the implicit time update 
in Sec. 13.31 


3.1. Advection in Space 


We employed a st andard high-resolut ion-shock-capturing scheme and utilize t he HLL (Harten- 
Lax-v an Leer) scheme ( Harten et alll983 ) to evaluate the numerical fluxes in space ( Kuroda fc Urneda 
201fll ) . As for the fastest and slowest characteri stic wave speeds of th e radiatio n field system (Eos. ((11) 


and ([5])), we again use the same definition as in lKuroda et al.1 ( 20121 ) (see also lShibata et al.1 ( 201ll )) 
and connect A ra d,thin and A ra d,thick smoothly via the variable Eddington factor x as 


Arad — n Arad,thin A r, A ra d,thick) 


(36) 


2 2 

where A ra d,thin and A ra d,thick are the wave speed in the optically thin and thick limit, respectively. 

To enforce the numerical flux of the radiation field in the opaque region asymptotically ap¬ 
proach to the diffusion limit, we evaluate the energy flux (E 9 U ) and the momentum flux (T^) 
as 


pO _ 

-^hii - 


A + F°-A_F°+eA_A + (Q° -Qj) 
A + — A_ 


and 


e 2 0 + q - A-Fj) + eA_A + (Qk - Qi) . 

Ull - - r - r -h (1 - c ) 

A_|_ — A_ 


2 ,n+F i R 


(37) 


(38) 


respectively (jAudit et al.ll2002l : lO’Connor &: Qttll2013l h Here, and ^l/r are ^ ie conservative 


variables and their corresponding fluxes, respectively, with L/R denoting the left/right states for 
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the Riemann problem. All the radiation-hydrodynamical variables are defined at the cell center. 
For those cell centered (primitive) variables, we use a monotonized central reconstr uction, w h ich has 
secon d order accuracy in space, and obtain left/right states at the cell surface (see iKuroda fe Umeda 
2010 , for more detailed explanation). After the reconstruction, all the characteristic wave speeds 
of the matter and radiation fields are evaluated. 

e is a modification parameter to fit the numerical flux to the diffusion limit, which we take as 

(39) 


where k is the total opacity and Ax is the grid width ([Audit et al 


2002 


O’Connor & OttJ 20131 ). 


3.2. Advection in Energy Space 


Regarding the advection term in energy space in all conservation equations (Eqns.(|4|), (|5j) 
and (1 16 II ). w e def ine the advection fluxes at the interface of the energy bin as the same as in 


Muller et ah (|2010h so that all energy integrated advection term will vanish. This can be achieved 
since both terms M a , appearing in energy and momentum conservation equations, and eq a ^ in 
number conservation equation are expressed in terms of linear combinations of radiation momenta, 
J, H a , L a P ... , and we therefore can define their cell surfaced values with an appropriate weighting 
function to suppress violation of all conservations simultaneously. 

For all orders of the radiation momenta X £ {J, H a , L a P the following conditions need to 
be satisfied for number conservation, 


ded £ X = 0 




Xi+l/2 — Xi_ 1/2 

A Ei 


= 0 , 


(40) 


and for energy and momentum conservations, 

J dsd £ (sX) = j de (ed £ X + X) = 0 = 


S Ae < 


_ x m /2 - 1/2 , v 

£«-7-T Aj 

A £i 


= 0, (41) 


respectively. The RHS of the above equations represent the finite difference expressions with i and 
i + 1/2 denoting z-th energy bin and the interface between z- and (z + l)-th energy bins, respectively. 
A £i is energy grid width Ae* = e i+ i/ 2 — £j_ i/ 2 - It is straightforwad to show that Ea. flOj) can be 
automatically satisfied for any cell-surfaced quantities as long as they vanish at the outer boundary 
in the energy space. By introducing a definition X i+1 / 2 = Xf + A^ 1 , the RHS of Ea. (|TTj) can be 
expressed as, 


^ [—(£»+! - Si)xl - (Si - £i-l)Xf + A £iXi] = 0. 


( 42 ) 
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As in IMiiller et al.l Q.2010), we can get the solution of Eq.(42) as 

A £i 


X L = 


x R = 


^•z+l &l 

A £i 


-x&, 


°i —1 

where £* is a weighting function and is expressed as 

fa 

Ji+l/2 




& 


fa I fa 

U- 1/2 ~ H+ 1/2 

In this study, we used a “Harmonic” interpolation (<j = 1) for the energy density as 


(43) 

(44) 

(45) 


ft 


+1/2 


E te) 

e? 


l-r; 


+ !/2 / E 


W) 

-3 

; i+l 


r i+l/2 


(46) 


with r i+1 f 2 = (^i+i /2 — £i)/(£i+i — £*) (see IMiiller et al. (120101 ) for more details). Like these, just 
only by evaluating an appropriate radiation momenta / 2 at the energy bin surface, we can 
simultaneously vanish Eons. (1401) - (1411) numerically and maintain energy, momentum and number 
conservations. 


3.3. Implicit Time Update 

After the explicit update, we solve the following simultaneous equation (e.g., Eq. ([35]) ): 


,n+U _ Q(P' t+1 ) - Q* 


f(pn+I) _ 


At 


+ S 


adv.e V 


yn+l\ 


+ S um (P n+1 ) = o, 


(47) 


where Q, S a d v ,e and S„ m are expressed in terms of the primitive variables P 


P 

Ui 

s 

Y P 


E, 


0+) 
\v,z)i - 


(48) 


at (n + l)-th time step. Obviously, the baryon number density does not change in this step, i.e., 
p* = p n+1 . To get the solutions of the above simultaneous equation, we employ the Newton- 
Raphson iteration scheme with the inversion of the following matrix until a sufficient convergence 
is achieved; 

<9f(P J ) 


<9P 


OP J = -f(P J 


>i+i 


= P 1 + 8P 1 , 


(49) 

(50) 
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for I = 0,1,2., with the initial condition P° = P*. As for the convergence criteria for the 

Newton-Raphson iteration, we monitor 


ISP 1 

|P 7 | 


< tol. 


(51) 


where tol represents a tolerance and we typically set tol = 10 . 

We note that, even if we evaluate the advection terms in energy space explicitly in time as 
S a dv,e n 5 Eas. (l47l) - ([5TT) do not change except the term S a d v ,e(P ?1+1 ) in Ea. (l47l) is now moved to the 
explicit part. 

Our method for time update from n-th to (n + l)-th time step is summarized in order as follows 
and schematically drawn in Figure [H 


1. Geometrical Update. We first evolve all the BSSN and the gauge variables G={ 7 ij, Aij, cf>, 
K , r*, a, /?*} from n-th to (n + l)-th time step. 

2. Explicit Update. In the second step, all the advection in space (S a d v ,s) and gravitational 
source (S grv ) terms are added to obtain the fractional time-step values Q* and their consis¬ 
tent primitive variables P(G n+1 , Q* ) = (p,Ui,s,Y e ,E^ £) ,F { ^ £) .). Advection of total lepton 
number (Ea. ([20l) ) is also performed simultaneously to evaluate Y*. 

3. Implicit Update. Finally, quantities in the fractional time-step (Q*) are implicitly updated 
to those in the (n-(-l)-th time-step (Q ra+1 ) by the Newton-Raphson iteration until a sufficient 
convergence is obtained with a constraint that the local lepton fraction does not change (i.e. 
Y] n+1 = Y* since p n+1 = p*) in the diffusion region. 


Here we shortly summarize some technical details to achieve a sufficient lepton number conser¬ 
vation in practical core-collapse simulations. As we have already mentioned in Sec. 12.21 neutrino 
number conservation is formally satisfied by solving the energy-momentum conservation equations 
([l])-©. Especially, in the neutrino trapping regime, solving the energy conservation equation © is 
practically identical to solving the neutrino number conservation equation (1161) which was proven in 
Sec. 12.2.11 However, in the finite difference method, it does not guarantee a perfect match between 
equations ©-© and equation (|16l) due to discretization error. To minimize the difference, we add 
the following constraint to Eg. (1511) 

v/ _ w*| 

1 Y I 1 < to1 ’ ( 52 ) 

as another criterion to exit the Newton-Raphson iteration. Because of this additional criterion, the 
lepton number conservation is also satisfied when the Newton-Raphson iteration converges. Note 
that the local baryon number do not change in Eq. (l47l) . i.e., Y* = YJ n+1 should be met if the lepton 
number is conserved. 
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G n , Q" fc P»(G ra ,Q ra ) 

^ Geometrical updating part 


G ra+I ,Q" fe P ra (G n+1 ,Q 


Explicit updating part, Eq.(32) 


i n+1 


, Q* & P*(G n+1 , Q*) 


Implicit updating part, Eq.(33) 


Perform following matrix inversion iteratively 

<9f(P 7 ) 


5P 1 


dP 1 


= -f(P') 


I 


* m<toi & 


Y/ 


1 ! p/ , 

i set Q n+1 = Q 1 , P n+1 = P 7 


< tol 




171+1 /^7l+l 


, Q n+1 & P n+1 (G n+1 , Q" +1 ) 


Fig. 1.— A flowchart to visualize how to update variables (G, Q, and P) from n-th to (n + l)-th 
time step (see text). 


We also adopt another prescription in which K n+1 is used to achieve a better convergence for 
the Newton-Raphson method. In some cases, Yj happens to go beyond the range of the employed 
EOS table during the iteration, especially when the Jacobian matrix is not well evaluated. In those 
cases, we use F) n+1 to reset Yj as below 


Y J = Y , n+1 - K + Yl- 


(53) 


With this reset, we found that Yj hardly goes beyond the range of the EOS table and also does not 
take unreasonable value for the supernova core. In addition, the number of the Newton-Raphson 
iteration can sometime be reduced. This is because the reset value Yj automatically satisfies the 
lepton number conservation Y^ = Yj + Yj — Yj f = Y, n+1 which leads to the faster convergence. 

Without these two prescriptions, we observed that the central lepton fraction at bounce 
deviates maximally ~ 0.05 from the value immediately after neutrino trapping sets in ( p c > 
10 12 gem” 3 ). Note that these ad-hoc constraints are introduced just to find a more efficient path 
toward convergence in the implicit time update. The criterion for convergence is solely given by Eq. 
m • Since we do not include Eg. (1201) into Eg. (1471) . there is no inconsistency between the number 
of simultaneous equations f(P) and the unknown variables P n+1 . 
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4. Radiation Tests 

In this section, we show results of several simple test problems to check the validity of our Ml 
radiation transport code. Except for Sec.4.4, a flat spacetime is assumed throughout this section. 


4.1. Diffusion Wave Test 


We first perform a diffusion w ave test, by which we check the validity of our flux implementation 
in the diffusion limit. Following P ons et al.l (120001 ). a Dirac (5-function type radiation source is 
initially located at r = 0. Then we follow the diffusion of the source into the optically thick 
medium with zero absorptivity (k„ = 0) and high scattering opacity ( k s = 10 2 , 10 5 ). The source 
term in Eqs.([l])-([5]) thus becomes 


= -k s H». 


(54) 


The 3D computational domain is covered by 100 equidistant Cartesian zones in each direction 
(x € [—0.5,0.5]) for a model with k s = 10 2 . This corresponds to a Peclet number Pe = k s Ax = 1. 
While, in the other model with n s = 10 5 , we raise the nested level to 2 to achieve sufficient resolution 
at the center and to check that the nested structure does not interfere the radiation propagation in 
the diffusion limit. In this model, 32 3 numerical cells cover the base domain x € [—0.5, 0.5] with 2 
nested level. The minimum grid width and Peclet number at the center are thus Ax = 1/128 ~ 0.78 
and Pe = 10 5 /128 ~ 780, respectively. Since Pe is much higher than unity, e in Eg. (1391) approaches 
zero which enables us to check whether the radiation transport in the diffusion limit is solved 
appropriately. 


Th e analytic al solution of the zeroth and first order radiation momenta profiles at time T 
([Pons et al. 2000) is given as 


E(r,T) 

F r (r,T) 



(55) 

(56) 


respectively. From Figured it can be seen that our code can reproduce the analytical results quite 
well. 


4.2. Shadow Casting Test 


Next, we move on to a shadow casting test to check the ability of our code whether anisotropic 
radiation field can be ap propriatel y ev olved in the free streaming limit. The initial setup is es¬ 
sentially the same as in Kanno et al. (120131 ). We set a perfect absorbing region at x € [2,3] 
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Fig. 2.— Diffusion tests with n s = 10 2 (Pe = k s Ax = 1) in upper two panels and with k s = 10° 
(Pe > 780) in bottom twos. Our results (crosses) for the evolutions of the energy density (left 
panels) and radial energy flux (right panels) are compared with the analytical ones (solid lines). 
Note that the simulations start at T = 1 and 200 for models with k s = 10 2 and k s = 10 5 , 
respectively. The model with k s = 10 ;) is performed with two nested level structure and we plot the 
spatial profile of Pe in the bottom-left panel for reference. Note that we reduce number of plots 
for our results (crosses) to avoid messy overplots. 
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and { y,z } € [—2,2] inside the numerical domain x € [0,12] and {y,z} € [—4,4], We im¬ 
pose a constant radiation ( E,F X ) = (1,/max) from the left boundary. Numerical resolution is 
Ax(= Ay = A z) = 12/192 and we set / ma x = 0.999. Within the perfect absorbing region, the 
absorbing opacity is set as n a Ax ~ 10 10 , otherwise n a = 0. The source term in Eqs.dH)-© thus 
becomes 


S " = -K a ((j - j eq )u /* + H *), (57) 

with J eq = 0 and = (1,0,0,0). In Fig. [3] we show three different time snapshots (T=3.5, 7.5 



Fig. 3.— Shadow test at different time slices T =3.5, 7.5 and 15.5 from top to bottom. Color 
contour represents radiation energy density in logarithmic scale. The absorbing region is expressed 
by white-dashed line and the expected radiation shock front (i.e., x = T) is denoted by a vertical 
green line. 


and 15.5) of the radiation energy density in a logarithmic scale. The absorbing region is expressed 
by white-dashed line and the radiation shock front (x = T) is denoted by vertical green line. From 
this figure, we see that the absorbing condition works appropriately in our scheme and the radiation 
front propagates with the light speed v = 1. Furthermore, the region beyond the a bsorbi ng box 


i s bar ely contaminated by the radiation. These features are in good agreement with Kanno et al 

(hod i. 
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4.3. Propagation in Free Streaming Regime 


The next test is a spherical expansion of radiation held from a point-like source into optically 
thin medium. The aim of this test is to check whether our code using the Cartesian coord inates 
can maintain the sphericity of the held during the expansion. By following Just et al . (j2015l ). we 
define the radiation source S with radius r$ = 1.5 which centres at the origin of the 2D Cartesian 
domain with x € [—7.5, 7.5]. We also dehne the purely absorbing region A centered at x = (3.5, 0 ) 
with radius r „4 = 1. Inside those two regions, the absorptivity At a (x) and equilibrium energy density 
J eq (x) are set as 


«a(x) 


10 exp{ — (4a/ x 2 + y 2 /r s ) 2 } , x € S 
10 ,xe4, 


(58) 


and 


J eq (x) 


lir 1 ,xeS 
0 ,xgd, 


(59) 


respectively. In other regions, we set At a (x) = 0 and J eq (x) = 0. We neglect the scattering opacity 
k s (x). The source term in Eqs.dH)-® is thus the same as Eq. (l57|h Regarding the initial setup for 
the zeroth and hrst order momenta of radiation held, we assume an arbitrary dilute radiation held 
as 


E = 10 “ 9 
\F\/E = 10 _1 °. 


(60) 

(61) 


The 2D numerical domain is covered by 196 equally distant Cartesian zones in each direction with 
two nested levels, i.e. Ax varies from 15/784 to 15/196. 


In this test, we use following formula for the variable Eddington factor x (Levermore 1984) 


3 + 4F 2 

5 + 2V4 - 3 F 2 ' 


(62) 


Only in this propagation test, we take two possible options for evaluating the characteristic wave 
speeds of radiation held in the optically thin limit X±. Assuming a hat spacetime and using a 
following closure relation in the optically thin limit 


pfiu = E 


pv 

F h F k ’ 


the hrst one is ( Shibata et al. 201l|) 


(63) 


\S,i , ■ 

A ± = max/mm 


EF i 

Fi.F k 


±- 


F l 




( 64 ) 
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and the second one is (Skinner & Os trike r 2013) 


A| 0,i = \ ti l F± ^(4 - 3F 2 - ^4 - 3 F 2 ) + 2/T 2 (2 — F 2 — j /y/l - 3F 2 . 


(65) 


Here, = cos 6 is determined by the angle 6 which defines the orientation of the energy flux F 
relative to the interface normal x l . Note that, in the rest of sections, we use only Eg. (1641) for the 
free streaming part. 

Another thing to mention is that we limit the flux factor / = y/rfi FiFj /E to be less than a 
maximum allowed value / max by modifying the first order moment as 


Ft ->• min(/ max //, 1 )Ei, (66) 

in every time step. Note that setting / max < 1 means that we add contribution from the isotropic 
radiation pressure P t ? kick according to Eq.([ 6 ]). The reason of this modification will be discussed later 
in this section. 

In Fig. [U we show results for three different models taking different values for / max and using 

s /so 

different evaluation formulae for A ± ' . Upper two panels (a,b) and bottom one (c) are models using 

Eq. (jMD and Eq. (1651) . respectively. In panels (a) and (c), we take / max = 1, while / max = 0.93 is 
taken in panel (b). The color represents the energy density E in a logarithmic scale at four different 
time slices (T = 1,3,5,7). Inner two squares represent the boundary of nested structure and the 
dotted circle shows the purely absorbing region A. From models (a) and (c), which use / max = 1, 
we find non-isotropy which can be seen as corrugated patterns behind the radiation front, obviously 
with the exception of the absorbing circle A and its shadowing region. Furthermore in model (a) 
at T = 7, another remarkable violation of isotropy is seen at 45° away from the coordinate axises. 
Note that, from snapshot at T = 7 in panel (a), this violation seems to be associated with the nested 
structure. We, however, confirmed that this is not the artifact of that by performing a run with 
zero nested level. By comparing models (a) and (c), the violation seen at 45° is eliminated in model 
(c). Model (c) takes into account the propagation angle (//) relative to the interface normal more 
explicitly than model (a) when we evaluate A®°. We therefore consider a more accurate evaluation 
in the characteristic wave speeds \± for the numerical flux (Eqs (l37l) - (l38|) f is the key to follow the 
radiation in optically thin limit properly. While in panel (b), in which we use / max = 0.93, we do 
not see any spherical symmetry breaking. 

We also apply this test to a supernova core profile. We fix the hydrodynamical background, 
which has a typical core profile at T p b ~ 100 ms, and follow only the propagation of neutrino 
radiation. The neutrino transport part is identical to our practical calculation reported later in 
Sec. [5] and all relevant neutrino-matter interactions, gravitational redshift and Doppler shift terms 
are taken into account. There is thus a transition from optically thick to thin regime. We calculate 
two models with different values for / max . A ^' 1 (Eq l64j) is used in both models to evaluate the 
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Fig. 4.— Time snapshots of spherical explosion test in the optically thin medium with different 
sets of (/max, A-)-). Models (a) and (c) (top and bottom rows) are models with taking / max = 1, 
while /max = 0.93 is taken in model (b) (middle row). Eqs. (f64l) and (165)1 are used in models (a-b) 
and (c), respectively. Color scale represents radiation energy density on a linear scale in arbitrary 
unit. Inner two squares represent the boundary of nested structure and the dotted circle shows the 
purely absorbing region A. 
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characteristic wave speeds. In Fig. 0 we show the radial profiles of (electron-type) neutrino 
luminosity for two cases: one is calculated with / max = 0.999 (black diamonds) and the other is 
with / m ax = 0.93 (red filled triangles). As can be clearly seen, beyond the shock position (R > 100 



Fig. 5.— Arbitrary normalized surface integrated local neutrino luminosities measured by an 
Eulerian observer. Background profile is a stationary and spherically symmetric SN core taken at 
post bounce time 100nrs. Models are with two different maximum allowed flux factors, / max = 0.999 
(black diamonds) and = 0.93 (red filled triangles). 


km) where it is optically thin regime, a large scatter is seen at R > 300 — 400 km for model with 
/max = 0.999. Spatial location of those highly deviated radiation are again concentrated along 
~ 45° away from the coordinate axises as the same as in the previous simple propagation test 
(panel (a) in Fig. [4j. Meanwhile, in model with / max = 0.93, there is little deviation and the 
surface integrated local luminosity stays nearly constant in the transparent region as it should be. 


From these two tests, we are currently enforced to set a ceiling value for the flux factor / < 
/max = 0.93 in the practical core-collapse simulation to follow spherical-like propagation stably. 
Since / max = 1 is the physically correct value, one may suspect that our modification can be an 
obstacle for comparison of the emergent radiation profile with pr eviou s studies. From previous fully 
relativistic Boltzmann transport simulation (iLiebendorfer et ajj 20011. however, it was shown that 
the flux factor higher than, e.g., ~ 0.93 appears only beyond R > 500 — 1000 km. Therefore even 
if we evaluate the emergent neutrino profile at the same radius R = 500 km as previous studies, 
it is altered only slightly ~ a few % and the modification cannot be a significant obstacle in this 
study. Obviously, the SN hydrodynamics itself is not affected by our artificial treatment since it is 
applied only to the optically thin region. 






- 24 - 


4 . 4 . Gravitational Redshift and Doppler Shift 


To 


check the energy coupling terms for gravitational redshift and Doppler shi ft, y/jad e (eMuri n') 


and '/y adJeM^ j/,) in Eqs.flU)-©, we repeat the same tests performed in lMtiller et al.1 (j201ol ): 


O’Connor! (2014). We consider the propagation of radiation from a sphere with radius R = 10 km 


through a curved space-time and sharp velocity profile in spherical symmetry. Regarding the sharp 
velocity profile, we take the following one 


u r 


0, r < 70km 

< -0.2c ( - 7 k °r ) , 70km < r < 80km 
- 0 . 2 c (fc ) 2 , r > 80km. 


(67) 


Here, u r is the radial component of Uj. As for the curved space-time, we take the same energy density 
profile for matter which is used in the second test problem in Sec. 14.31 and solve the Hamiltonian 
constraint to obtain the conformal factor ip and lapse a. When we solve the Hamiltonian constraint, 
we assume zero velocity profile (c,; = 0 = fij) and the conformally flat approximation ( 7 \j = ip^'pij = 
As for the initial neutrino profile within the radiating sphere, we take zero chemical potential 
/T, = 0 and a temperature of 5 MeV, i.e., (e) ~ 15.7 MeV, in the free streaming limit E = y/ 7 ijF^Fi. 
Here (e) is the mean energy of neutrinos as measured in the comoving frame. Outside of the central 
radiating sphere, we choose an arbitral dilute radiation field with (e) ~ 15.7 MeV. We do not 
evolve neutrinos inside the radiating sphere and follow the propagation only at R > 10 km. The 
3D computational domain is a cubic box with 3000 km width (i.e. the outer boundary is at the 
radius of 1500 km from the origin) and nested boxes with 7 refinement levels are embedded. Each 
box contains 64 3 cells and the minimum grid size near the origin is Ax = 366m. As for the energy 
grid of the neutrino radiation field, we use logarithmically spaced 20 energy bins (N e = 20) which 
center from e = 1 MeV to 300 MeV. We calculate two models with and without the sharp velocity 
profile. Both models are calculated in the curved space-time. 


For given velocity and space-t ime profi les, the free stre aming neutrinos propagate with obeying 


the following relations (M uller et al 


201C; 


O ’Conno r 2014) 


l^a(l + u r )(e) = const, ( 68 ) 

y^yar 2 'y rr j deF r = L eu i = const, (69) 

where W is the Lorentz factor and v r = u r /uP. = ip 6 and 'y rr = are the determinant and 

rr-component of the three metric, respectively, in the conformally flat approximation. L eu i is the 
luminosity measured by an Eulerian observer. Eg. (1691) is derived from the stationary solution of 
Eq. d4]) with neglecting the source terms. 

In the top panel of Fig. [ 6 j we show the background profiles for ip, a and 1 — u r . (e) and L eu i 
are plotted in the middle and bottom panels, respectively. In the middle panel, our results (filled 
triangles) and analytical ones (solid lines) are plotted. The analytical expression for (e) is derived 
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Fig. 6.— In the top panel, we show background profiles for if;, a and 1 — u r . The mean energy 
(s) and Eulerian luminosity L eu i profiles are plotted in the middle and bottom panels, respectively. 
Two test cases, with (u r ^ 0) and without (u r = 0) the shock profile, in the curved space-time 
are plotted. In the middle panel, our results are denoted by filled triangles and the solid lines are 
analytical results. Note that (e) for the model with the shock profile is shifted upward with IMeV 
to avoid the overlapping of plots. 
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from Eg. (1681) where the constant in the right hand side is evaluated at R = 10 km. Numerical 
results are taken after they settle down in nearly stationary. 

As seen in the middle panel, our results reproduce the analytical ones quite well except at 
the shock front. We consider that the less agreement at the shock front originates from the coarse 
spatial resolution. The finite difference expression for the term 'S/ a up in the energy coupling terms 
cannot capture the sharp velocity profile with enough accuracy. From the bottom panel, we can 
find the Eulerian luminosity stays nearly constant with a few percent deviation beyond R > 20 km. 
Note that the zig-zag pattern seen in L eu \ is an artifact of the mesh refinement. 


5. Core Collapse of a 15 M & Star 


In order to confirm the validity of our new supernova code, it is of primary importance to make 
a detailed comparison with the previous published results. We employ the data from 1D-GR neu¬ 
trino transport si mulations (ILiebendorfer et alJl2005l : ISumiyoshi et al.ll2005l : iMiiller et al .1120101 ) and 
from 2D-GR ones ( Muller fc Jankal 20141 ). We chose these models because all of them took the same 


proge nitor and employed similar microphysics in the GR CCSN simulations. ILiebendorfer et al 


( 200511 presented detailed comparison of two independent numerical codes, AGIL E-BOLTZTRAN 
( Liebendorfer et al. 2004 ) and VERTEX-PROMETHEUS ( Ramon fe Janka 2002 ). Since their re¬ 
sults are available onlinqj, we are able to make a detailed comparison with their data set. AGILE- 
BOLTZTRAN solves the GR Boltzmann equation with the S n method in spherically symmetric 
Lagrangian mesh, whereas VERTEX is an Eulerian code that solves the moment equations of a 
model Boltzmann equation by the VEF method in the Newtonian hydrodynamics plus a modi¬ 
fied GR potential (VE RTEX-PROMETHEUS) and also in the c o nform ally-flat GR hydrodynamics 
(VERTEX-CoCoNuT, Dimmelmeier et al. (2002); Miiller et al. ( 2010 )). Our code is rather simi¬ 
lar to VERTEX-CoCoNuT than AGILE-BOLTZTRAN except for the different geometrical solvers 
and the different coordinate systems are adopted. In the following, we label t he resu l ts o f AGILE- 
BOLTZTRAN as “ABG15”, of VERTEX-PROMETHEUS as “VXG15” and of lMiiller et J (I2OI0I ) 
as “BMG15’fl. 


5.1. Numerical Setups 


We employ a 15M 0 progenitor (model “sl5s7b2’ : 
collapse, bounce and initial postbounce phase up to 


Wooslev & Weaver 19951) and follow core 


~ 50 ms. We use the EOS of Lattimer & 
Swesty (1991) (LS EOS) with an incompressibility parameter of K = 2 20 MeV . Even t hough our 


choice of K = 220 MeV is different from the one (K = 180 MeV) used in 


Liebendorfer et al 


(20051); 


5 We used data which can be downloaded at http://iopscience.iop.Org/0004-637X/620/2/840/fulltext/datafiles.tar.gz 
6 We read their data from the figures digitally and plot them in this paper. 
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Sumivoshi et al.1 (1200.1 ) ; iMiiller et al.1 (|20ld ): lMuller fe Jankal (120141 ) . rThonipsoii et al.l (J2003j) showed 
that differences in the neutrino profiles among models with the different K of LS EOS are a few % 
around core bounce and less than ~ 10 % for the first ~ 200 ms after bounce. We thus consider 
that the different choice of K barely disturbs the aim of our comparison study. 


_ As for neutrino opacities, the standard w eak interaction set inlBruennl (119851 ) and lRampp fc .Tanka 

( 2002 ) plus nucleon-nucleon bremsstrahlung ( Hannestad fe Raffeltl 199^ ) is taken into account (see 
Table 1 and Appendix [A] for more detail). For simplicity, we neglect higher harmonic angular de¬ 
pendence of the reaction angle when we calculate the source terms for neutrino electron scattering, 
thermal pair production and annihilation of neutrinos, and nucleon-nucleon bremsstrahlung (i.e., 
^nes/tp and C \e), nes/tp are Set t0 be °> See A PP S - D and [AT]). 


Process 
nv e e~p 
pv e e + n 
v e A o e~ A' 
iyp -H- vp 
vn <->• vn 
vA <H- vA 
-H- ve± 
e~ e + VV 

NN ++ vvNN 


reference 

Bruenn 11985 ); Rampp fc Janka 12002) 
Bruenn 11985) ; Rampp &; Janka 12002) 
Bruenn 11985) ; Rampp &; Janka 12002) 
Bruenn 11985) ; Rampp &; Janka 12002) 
Bruenn 11985) ; Rampp &; Janka 12002) 
Bruenn 11985) ; Rampp &; Janka 12002) 
Bruenn 11985 ) 

Bruenn 11985 ) 

Hannestad Sz Raffelt (1998) 


summarised in 
AppOO 
AppjAT] 
AppjAT] 
AppjA^l 
AppJXl 
AppD 
AppjA3] 
ApplAAl 
AppJAA] 


Table 1: The opacity set included in this study and their references. Note that v, in neutral current 
reactions, represents all species of neutrinos (v e , V ei v x ) with v x representing heavy-lepton neutrinos 
(i.e. Vn,v T and their anti-particles). 


In this study we investigate two models, 1DG15 and 3DG15, both of which are computed in 
the 3D Cartesian coordinates. While model 3DG15 is calculated without any artificial constraints, 
model 1DG15 is considered to mimic ID model by artificially suppressing the non-radial component 
of the flow velocity Ui as 


x J u 


Ui = 


V. 


(70) 


Although this artificial elimination coul d potent ially lead to the shift of the kinetic energy into the 
thermal one, our previous study (Kuroda et ah 2012) showed that the violation of the momentum 
constraint is negligible during the early postbounce phase (T p b < 100 ms where T p b denotes the 
postbounce time). 


The 3D computational domain is a cubic box with 8000 km width (i.e. the outer boundary is at 
the radius of 4000 km from the origin) and nested boxes with 5 refinement levels, at the beginning 
of calculation, to 8 refinement levels, when the central rest mass density reaches 5 x 10 13 g cm -3 , 
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are embedded without any spatial symmetry. Each box contains 64 3 cells and the minimum grid 
size near the origin at bounce is Ax = 488m. In the vicinity of the stalled shock front R ~ 100 km, 
our resolution achieves Ax ~ 3.6 or 7.2 km, i.e., effective angular resolution becomes 3.6km/100km 
~ 2° or ~ 4°, which is considered to be rather too coarse resolutio n to follow a nearly spherical 


structure in the Cartesian grids. Compared to recent 3D-GR studies (jOtt et al. 20121 ; IKuroda et al 


20141 ). the resolution is still approximately two times coarser. The time step size is At = 10 7 s 


which corresponds to ~0.06 in the CFL number. The main reason for taking such a small timestep 
is to get a rapid convergence during the iteration in our implicit time update part (see, Sec l3.3ll . 
which is not restricted by the CFL condition (of the explicit update in a pure hydrodynamics case). 
As for the energy grid of the neutrino radiation field, we use logarithmically spaced 20 energy bins 
(N e = 20) which center from e = 1 MeV to 300 MeV. 


5.2. Results 

In this section we start to make a detailed comparison firstly before bounce (section 15.2.11) 
and then after bounce (section I5.2.2|) between our code (models 1DG15 and 3DG15), AGILE- 
BOLTZTRAN (model ABG15), VERTEX-PROMETHEUS (model VXG15), and (partly) from 
VERTEX-COCONUT (model BMG15). In the following, we call the results from the latter three 
codes as the reference results. 


5.2.1. Before bounce 


In the upper panel of Fig. [7J we plot the central (matter) entropy s, electron fraction Y e and 
the total lepton fraction Y) = Y e + Y u as a function of the central (rest-mass) density p c for model 
1DG15 (black), ABG15 (blue), and VXG15 (green), respectively. The lower panel of Fig. 1 shows 
the evolution of p c and the deviation of the ADM mass AMadm from its initial value. Here Madm 
is given by 


Madm = 


+ 



(<Sb + J deE^^je- 
e- 4 >[s l + J deF^j 


p 5(f> / 9 \ 

+ — ( AnA* - -K 2 - fiR^e-M) 
16vrV 13 3 ' J J 



(71) 


where the second line denotes energy loss due to momentum and neutrino energy flux through the 
numerical boundary and n represents a unit normal vector to the surface element da. In the above 
surface integration, we ne glect energy loss due to g ravitational wave emission since it is negligibly 


small 10 11 M«c 2 , e.g., 

Scheidegger et al. 


201C)) compared to the violation of the ADM mass 

in the CCSN environment (e.g., 

Kotake 

(2013 

) for a review). 


The top panel of Fig. 1 shows that the neutrino trapping starts p c ~ 2 x 10 12 g cm 3 for 
model 1DG15 and the lepton fraction remains constant with Y) = 0.323 afterward. The evolution 
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Fig. 7.— Upper: the central (matter) entropy s, electron fraction Y e and total lepton fraction 
Yi = Y e + Y v as a function of central density p c for model 1DG15 (black), ABG15 (blue) and 
VXG15 (green). Lower: Comparison of postbounce evolution of the central density (solid) between 
the four models and the deviation of the ADM mass AMadm (dash-dotted) for our two models 
with (model 1DG15) or without the artificial elimination of the non-radial velocity (model 3DG15, 
see text). 
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of our c entral Y P and Yi (bl a ck lines) is quantita tively in good agreement with VXG15 (green line, 
see also lMiiller et al.1 ( 2010 1: iBuras et al.1 ( 2006bl ll and ABG15 (blue line). As already explained in 
Sec. 12.21 we solve the advection equation (Eg. (l20l) l of the total lepton number (density) pY/ as a 
constraint to ensure the lepton number conservation. Because of the treatment, the evolution of the 
lepton/ electr on fraction is in excellent agreement with the reference results. As already pointed out 
by O’Connor (20141), this also relies on the accurate implementation of inelastic neutrino-electron 
scattering, energy-bin coupling, and the appropriate closure relation. 

After the core deleptonization ceases (i.e., Yj stays nearly constant with increasing central 
density), the inner core evolves almost adiabatically and the entropy remains nearly constant as it 
should be. A small breaking of the adi abatici ty (d ecrea se by 3.8% in the central entropy) is seen in 
model 1DG15 before bounce (see also O’Connor (2014)). However, the (time-a veraged) value s 


1.3 kp , baryon 1 is rough l y in g o od agreemen t with t he ref erence results ( see al so lLiebendorfer et al 


( 20051 ): ISumivoshi et al.l ( 20051 )) : IBuras et al.1 ( 2006bl h and lMiiller et al.l ( 20100 ) and this would not 
have a big impact on the subsequent core evolution due to the short simulation time in this study. 

As for the total energy conservation (bottom panel of Fig.3), it is maximally violated with the 
amount of AMadm 4 x 10 4 Mq for model 1DG15 (i.e. r*j 8 x 10 50 erg). The v iolation at bounce is 
slightly worse than that (Aadm ~ 5x 10 50 erg) of the VERTEX-CoCoNut code ( Muller et al. 2010 ). 
which should be improved and kept much smaller in more precise CCSN modelling. As one would 
anticipate, the violation is bigger for model 1DG15 (dashed black line) with the artificial elimination 
of the non-radial velocity than for model 3DG15 without (dashed red line). In our previous study 
(with the gray Ml scheme, Kuroda et aL, (2012)). the violation of the ADM mass is typically one 
order-of-magnitude smaller than that for the corresponding 3D model with (approximately) twice 
higher resolution. Because the computational time for the implicit update in the current code is 
very expensive (using our best available resources), we are now forced to employ a quite coarse 
resolution. It is important to clarify how the violation of the total energy conservation would be 
improved with increasing numerical resolution. We leave this for future work. 
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10 
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Fig. 8.— Neutrino distribution function /(z/, e) (filled triangles) and Fermi-Dirac distribution func¬ 
tion at equilibrium (solid lines) at the innermost mesh for model 1DGR. Lines and triangles are color 
coded according to the infall phase. Note that f(u, e) for anti-electron neutrino (P e ) is multiplied 
by 10° for comparison. 
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Fig. [8] shows a spectral shape of the neutrino distribution function /(V, s) (filled triangles), 

J (u,s) 


/(^e) = 


47T£ 3 


(72) 


which is estimated at the innermost grid point of model 1DG15 when the central density p c reaches 
lO 11,12,13 ’ 14 g cm" 3 , and at bounce, respectively. Solid curves represent the Fermi-Dirac distribution 
at equilibrium. From the left panel, it can be seen that /3-equilibrium is achieved for electron-type 
neutrino (zz e ) at 10 12 gcm~ 3 < p r < 10 13 gcm~ 3 , which is consistent with the neutrino trapping 
density as shown in Fig. [TJ In Bruenn (1 9851 ), the /3-equilibrium for v e was obtained after p c = 
2.46 x 10 12 g cm -3 in their corresponding model to ours (“standard” model). Note that at bounce, 
electron type neutrinos, in the low energy range < 5 MeV, slightly violate the ferrni blocking, 
i.e., the distribution function /(zz e ) exceeds one. The excess, however, is < 10 -5 and we consider 
that it is a negligible amount. As shown from the middle (P e ) and right panel ( ux ) of Fig.4, other 
neutrino specie s are thermalized only af t er p r e x ceeds 10 14 g cm~ 3 . Thes e features are quantitatively 
consistent with iBruenn fc Mezzacappal (|1997l ) ; iRampp fc .Tankal (120021 ) . It may be surprising that 
regardless of the use of different EOS and different hydrodynamics codes, the trapping density 
of v e (ptrap = 2 ~ 3 x 10 12 g cm" 3 ) in modern sim ulations (e.g. , top panel of Fig.3) is in good 
agreement with the pioneering work in the 1980’s (Bruenn 19851 ). It should be also noted that 


rates (Laneanke & Martinez-Pinedo 

2000; 

Juodagalvis et al. 

2010) need to be implemented as in 

Langanke et al. ( 

2003 

); 

Lentz et al. ( 

2012a 

b). 


So far, we have shown that our Ml scheme can capture several important phenomena regard¬ 
ing deleptonization, such as the neutrino trapping and the conservation of the lepton fraction in 
the diffusion region. As already denoted in Sec. 12.11 this is not trivial in the finite difference 
method especially when one transports conservative radiation moments (E^, F^.) instead of the 
corresponding comoving variables of (77( £ ), 73( e ).). A key is to find an appropriate Eddington factor 
Xt £ ) by which the neutrino energy flux approaches F/ e \ l —>• 4/3 E^v 1 in the diffusion limit (see, 
Ea. (l26l) h Since this relation can be achieved only when = E^yi /3 holds, X( e ) and F should 

approach 1/3 and 0 (e.g., our closure relation (Eq. (|6j))) in the limit, respectively. 

To show that both F = /J 2 ^ and F^ 1 actually approach 0 and 4/3 E^v 1 in 

the opaque region, we plot in Fig. [9] the radial profiles of (F r )/(E) (solid lines) and (H r )/(J) 
(dash-dotted lines) at different p c . Here (X) = f deX represents the energy integration of X. 
At p c = 10 12 g cm -3 , both solid and dash-dotted green lines almost coincide. However, as the 
infalling matter velocity comes closer to be relativistic ( p c > 10 13 g cm" 3 ), both lines start to split 
especially within R < 70 km. In the central region ( R < 70 km), (H r )/(J) approaches 0 towards 
the center, whereas {F r )/(E) becomes negatively large with the peak being around R ~ 30 km and 
then converges to zero to the center. We also plot the radial velocity of matter ( V r ) measured in 
the Eulerian frame and multiplied by 4/3 (blue diamonds) at p c = 10 14 g cm" 3 . 


Because of our appropriate evaluation for the Eddington factor, the flux factor measured in 
the Eulerian frame (approximated here by ( F r )/(E)) nicely matches with 4V//3 in the optically 
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Fig. 9.— Radial profiles of (F r )/(E) (solid lines) and (H r )/(J) (dash-dotted lines) for electron-type 
neutrino (u e ) when the central density reaches p c = 10 12,13,14 g cm" 3 in model 1DG15. V r denotes 
the radial velocity of matter (see text). 


thick region. This neutrino advection is essentially important for the radiation energy (Er e \) to 
move with the same velocity v l with matter in the opaque region (see, Eqs. (I26I) - (I29I) 1. He re w e 
shortly comment on the definition of the flux factor F. In our previous study (IKurod a et al.ll2012 ~) 


we em ployed the definition F = yJ^FiFj/E which is one of the candidates for F (jShibata et al. 


201 11. As can be clearly seen from the split between the solid and dashed-dotted lines in Fig. [9j we 


show that our previous choice of the flux factor is not adequate because the optically thick medium 
moves (albeit mildly) relativistically in the collapsing core. 


5.2.2. After bounce 


In Fig. m we show the radial profiles of various quantities (the rest mass density p, radial 
velocity v T , matter temperature T, entropy s, averaged neutrino energy E v and electron/neutrino 
fraction Y e jY\) at the selected time slices shortly after bounce. Here, we define the average neutrino 
energy E v as 


_ f de e 3 / (z/, e) 
~~ fdee 2 f{y,e )' 


(73) 


When the central density exceeds nuclear saturation densities (top left panel of Fig JTOl) . the 
core bounces because of the strong repulsive force of nuclear matter. Then the bounce shock 
is formed (green line, left middle panel of Fig llOl) which propagates outward with dissociating 
infalling heavy nuclei into free protons and neutrons. Production of enormous free protons/neutrons 
significantly enhances electron capture process, e~p —>■ v e n, behind the stalling shock. Immediately 
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Fig. 10.— Radial profiles of the rest mass density p , the radial velocity v r normalised by the speed 
of light c, the matter temperature T, the entropy s, the averaged neutrino energy E u and the 
electron/neutrino fractions Y e /Yi at selected time slices shortly after bounce T p b = 0, 1, 3, 5 ms. 
Note in this plot that profiles only along the x-axis of model 1DG15 are shown (because model 
3DG15 shows almost the same profiles). 
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after bounce, those high-energy neutrinos (top right panel) are still trapped inside the optically 
thick medium behind the shock. The medium, however, quickly becomes transparent to neutrinos 
and those neutrinos are liberated suddenly as a burst. This neutronization burst enhances further 
electron capture behind the shock due to the continuous deviation from the /3-equilibrium, leading 
to a characteristic trough in the Y e profile seen at R ~ 50 km behind the shock (right bottom 
panel of Fig llOD . Due to the energy loss by the photodissociation of the iron nuclei and the rapid 
neutrino leakage, the bounce shock stalls at R ~ 70 km within T p b ~ 3 ms (left panels of Fig llOD . 


Such dynamical f eatures are commonly seen in the reference models (jLiebendorfer et al 


2005; 


M uller et al. 20101 ). This indicates that our Ml scheme can capture the basic behaviours of the 


neutrino propagation from the optically thick to thin medium (otherwise it would result in either 
the absence or the different position of the Y e trough). 



Fig. 11.— Radial profiles of the (electron-type) neutrino energy flux (solid lines) and the radial 
velocity (dash-dotted lines). Numbers beside each line (1,2,3..) correspond to time slices, which 
are denoted in the lower right part with T p b in ms (0.7, 1.7, 2.7..). For the radial velocity in left 
panel, we plot only the first four time slices. 


How the neutronization burst is produced is more clearly depicted in Fig. [TT] where we plot 
the radial profiles of the neutrino ( v e ) energy flux (solid lines) and the radial velocity of matter 
(dash-dotted lines). At T p b = 0.7 ms, enormous neutrinos are still trapped and confined behind the 
shock, which is shown as a sharp peak in the energy flux around 10 < R < 20 km. At T p b = 1.7 ms, 
these neutrinos overtake the shock front because of the lowering opacity outward. Then the pulse 
of the neutrino burst propagates freely to the optically thin region (time label larger than 4) and 
eventually emerge s out of the computationa l domain (labels 9 and 10 ). Th ese profiles are consisten t 
with the results of lBruenn fc Haxtonl ()199ll ). iThompson et al.1 (120031 ) . and lR.ampp fc Jankal (120021 ') . 


Now we move on to make the code comparison in the postbounce phase. Figll2l compares the 
evolution of the shock radii for five models; two from our code (1DG15 (black line) and 3DG15 
(red line)), one from AGILE-BOLTZTRAN (ABG15, blue line), and two from VERTEX with 
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O 



Fig. 12.— Time evolution of the shock radius for five different models 1DG15, 3DG15, ABG15, 
VXG15, and BMG15. As for our models, we plot the angle averaged values. 


PROMETHEUS (VXG15, green line) and with CoCoNuT (BMG15, light blue line). Note that the 
final simulation time (T s j m ) is 50 ms for model 1DG15 (black line), whereas T s i m is 32 ms for model 
3DG15 (red line) simply limited by our available computational resources. 


The deviation seen in model 3DG15 (red line in Fig. fT2l) from the rest of the ID models is 
remarkable especially after T p b ~ 5 ms. This is because the bounce shock expands more energeti¬ 
cally in 3D pus hed p rimarily by prompt convection behind the shock. Using the same progenitor, 


Muller et al.l (l201 2bl) showed that the average shock radius becomes larger in 2D (their model G15) 
than in ID (their model G15-1D) at T p b ~ 70 ms because of the hot-bubble convection starts 
which is seeded during the deceleration of the prompt shock. The Cartesian coordinates have the 
intrinsic quadrupole perturbation and it affects significantly on the growth of the prompt convec¬ 
tion. The post bounce time, when the significant multi dimensionality appears, thus differ between 
our model and the ones using the spherical polar coordinates. In addition, our coarse numerical 
resolution might also lead to the earlier appearance of the initial convection because of even larger 
seed perturbations. 


The larger shock radius i n model 3DG1 5 than that in 1D G15 is a lso co n sistent with our prev ious 
result with leakage scheme (IKuroda et al.1 (|2012l ) , see also ICouchl (120131 ) ; lHanke et al.l (120121 ) for 
extensive discussion about the dimensional dependence on the postbounce dynamics). Now let 
us focus on our (pseudo-)ID model (black line in Fig. fl2l) . The shock radius of our code is in 
good agreement with the reference results exceptionally before T p b < 20 ms, whereas the shock 
radius tends to be smaller until the end of the simulation time. We consider that the difference 
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could primarily come from the use of the Cartesian coordinates with low numerical resolution and 
not from the neutrino-matter interaction terms. This is because our calculation in ID spherical 
coordinates with using the same neutrino-matter interaction terms shows a good agreement in the 
shock evolution (see, Appendix [B]). We give a more detailed discussion elsewhere below. 



Fig. 13.— In the clockwise direction from the top left panel, we show radial profiles of the (angle 
averaged) rest mass density p, entropy s, electron fraction Y e , and radial velocity V r at T p b = 3 ms 
for models “(1DG15)”, “ABG15”, and “VXG15”. 


In Fig. [13] and [141 we show various quantities from our code (labeled by “(1DG15)” and 
“(3DG15)” only in Fig. ED, AGILE-BOLTZTRAN (“ABG15”), and VERTEX-PROMETHEUS 
(“VXG15”) at two different postbounce times. 


At 3 ms after bounce fFig |13l) . the (angle-averaged) radial position of the stalled shock (bottom 
left panel) is R ~ 70 km for model 1DG15 (thick solid line). As seen, the velocity profile matches 
more closely the profile from AGILE-BOLTZTRAN (ABG15) than from VERTEX-PROMETHEUS 
(VXG15). This is also the case for the profiles of the density (top left panel), entropy (top right 
panel), and Y e (bottom right panel). It should be noted that the more recent res ults from the 
VERTEX-PROMETHEUS code with an improved GR potential ( Marek et ah 2006 ) agree very 
well with the AGILE-BOLTZTRAN code, hence with our code. Therefore we mainly compare to 
model ABG15 in the following. 


Looking at Fig. [13] more closely, one can see that the profiles of our entropy (top right panel) 
and Y e (bottom right panel) differ appreciably from model ABG15 especially in the region behind 
the stalled shock (R < 70 km) and above the unshocked inner core (R > 10 km). Let us remark that 
the early postbounce evolution starting from the shock formation, followed by the emergence of the 
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neutronization burst, until the shock stall is numerically most challenging. The code difference from 
the shock capturing scheme as well as the treatment of GR, the accuracy of the neutrino transport 
schemes could potentially impact the radiation-hydrodynamics evolution at the transient phase 
(i.e., 3-5 ms after bounce). Another remarkable difference is seen in the electron fraction. Among 
three results, only our result shows negative bump in Y e profile at 10 < R < 20 km. However this 
bump disappears in our ID spherically symmetric test problem in which we solve the advection 
terms in energy space S a d v ,e explicitly in time (see, Appendix [B] for more detail). 

Regarding the shock capturing, AGILE-BOLTZTRAN uses an artificial viscosity type with 
the second order accuracy in space, whereas our code employs the approximate Riemann solver 
(HLLE scheme like the VERTEX code) with the second-order accuracy in space both for radiation- 
hydrodynamical and geometrical variables. Thus the hydrodynamics part of our code is slightly 
more accurate than AGILE-BOLTZTRAN. On the other hand, the use of the approximate closure 
relation apparently falls behind the Boltzmann code especially in the semi-transparent region. 
Above all, the use of the Cartesian coordinates, which is very common in full-GR simulation^ makes 
the comparison to the genuine “ID” results of the reference models (based on the ID Lagrangian 
code (AGILE) and the multi-D code using the polar coordinates (VERTEX)) even more challenging. 



Fig. 14.— Same as Fig. [13l but at T p b = 50 ms for models “(1DG15)”, “ABG15” and “VXG15”. 
“(3DG15)” at T p b=32 ms is also plotted as a reference by a thin solid line. 

At Tpb = 50 ms (Fig {H]> , the differences between our pseudo-ID model (1DG15, thick line) 


See, however, iSanchis-Gual ct al . (120141 1 about recent report of the code development of numerical relativity in 


the polar coordinates. 
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and the reference results still remain to be seen rather remarkably in the postshock region (R > 100 
km), however this is not surprising given the different shock evolution (Fig. fl2l) . Here we consider 
that the numerical resolution in the postshock region sensitively affects the shock evolution. In the 
current resolution, the typical grid size of our nested box is ~ 7.8 km at 120 < R < 240 km (~ 4° 
resolution). As shown in Figll2l (red line), the deviation of the shock radius from the reference 
models becomes remarkable at T p b > 20 ms, which roughly coincides with the time when the shock 
reaches to the coarser level of the nested grid. There the shock front is resolved only by a few grid 
cells. We consider that at least a factor of two or more higher resolution is required to reproduce 
ID results, i.e., to recover the sphericity of the system in the Cartesian coordinates. But this is not 
an easy task as we will discuss the code performance in section [6j 


From Fig. El one can also see that the profile of the density, velocity, entropy, and Y e at 
the central region (R < 60 km) agrees with the reference results roughly within an accuracy of 
~ 10%. Given the insufficient numerical resolution, this good agreement in the semi-transparent 
region woul d su p port the validity of th e pres cribed closure scheme, which is in line with the recent 
results by O’Connor (20141) and Just et al. (2015). Fig El also shows that the profile of model 
3DG15 (thin solid line) differs from that of model 1DG15 (thick solid line). 


The entropy profile for model 3DG15 (thin line) behind the shock (R < 130 km) becomes 
more flat compared to the ID models, which is due to convective mixing behind the shock. As a 
result, the profiles of Y e and entropy in model 3DG15 become slightly closer to the reference results, 
though the similarity is just a coincidence and is not meaningful. 


FigEl shows the comparison of the neutrino luminosity L u which is the surface integral of the 
energy integrated comoving energy flux . through the surface of cubic box with 1000 km width 
(i.e. approximately 500 km from the origin). The peak luminosity, L Ue : ,peak = 4.1 X 10 53 erg s *, 
well agrees with 3.85 x 10 53 (ABG15), 3.8 x 10 53 (VXG15), and 4.3 x 10 53 (BMG15), respectively. 
After the neutronization phase when the ID core enters to a quasi-static phase (T p b > 20 ms), 
the anti-electron and heavy-lepton neutrino luminosities show quite consistent behaviors with the 
reference models and the differences between our results and, e.g., ABG15 are as small as ~ 10 %. 


Regarding the v e luminosity, we see a systematically lower value than the other three reference 
models. The difference reaches ~ 30 % (or ~ 10 52 erg s i) compared to ABG15 at T p b — 50 ms. 
As a consequence, u e and v e luminosities become comparable at T p b ~ 45 ms which is signi ficantly 
earlier than the reference values T p b ~ 70 ms in ABG15, BMG15 and also in a model of O’Connor 
(120141 ). This means that the lepton number loss from PNS is less than those previous studies since 
it can be measured roughly by the time integration of L Ve /e Ve — Lp e /ep e . However, our ID spherical 
test, albeit in the Newtonian limit, shown in Appendix [Bl does not exhibit such inconsistency 
and show quite reasonable neutrino profiles with a previous study. We therefore consider that the 
reason for less agreement, especially seen in L Ue , mainly comes from the spatial advection term in 
the Cartesian coordinates and not from the local neutrino matter interaction terms. 


We also compare the neutrino spectral difference in Fig. [T6l In the figure, we show time evolu- 
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0 20 40 

T P b ( ms ) 


Fig. 15.— Same as Fig. [12]but for the neutrino luminosity L v of electron type neutrinos (top panel) 
and anti-electron type and heavy-type neutrinos (bottom panel) for the five different models. 
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tions of the root-mean-square (RMS) energies of emergent neutrinos E Vt rm s mea sured at th e su r face 


of cubic box with 1000 km width. For E U:Tms , we use the same definition as in Liebendorfer et al. 
(2005) and it is defined as below. 


Ei/, rms — 


2 ' 


(74) 


fdef(v, e)s 4 

- V f def(u 

Again in Fig. [T6l we plot averaged value (E u rrns ) over the surface of the cubic box. From the figure, 



T pb (ms) 


Fig. 16.— Time evolution of RMS energies of emergent neutrinos E u rms for five different models. 
Our results are measured at the surface of cubic box with 1000 km width. 

we find good agreement with the reference codes and differences are less than ~1 MeV for u e and 
v e and <2 MeV for vx after the neutralization phase ceases T p b > 20 ms. In our model 1DG15, 
we see a spurious second peak in E VXirms profile around T p b ~ 13 ms. However, the second peak 
disappears in our model 3DG15 and we thus consider that it is most likely due to our artificial 
treatment for the matter velocity and is insignificant. 


6. Summary 

In this paper, we have presented our newly developed multi-D full GR neutrino-radiation- 
hydrodynamic code. The code was designed to evolve the Einstein field equation together with the 
GR radiation hydrodynamic equations in a self-consistent manner while satisfying the Hamiltonian 
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and momentum constraints. Using an Ml closure scheme, we solved spectral neutrino transport of 
the radiation energy and momentum in conservative way based on a truncated moment formalism. 
Beside the energy and momentum conservation, we paid particular attention to the lepton number 
conservation, especially in the neutrino trapping regime. 

We explained formally that the neutrino number is transported appropriately by solving only 
the energy and momentum conservation equations of the neutrino radiation field, especially focusing 
on the neutrino trapping regime. In addition, we showed that the advection terms in energy space 
are essential for reproducing the neutrino trapping. 

To validate our new numerical code, we first made bottom-line tests such as the diffusion test, 
shadow casting test and spherical propagation in free streaming regime to check the advection 
term in space. In addition, as for the important factor in CCSN simulation, gravitational redshift 
and Doppler shift are tested in a curved space-time with a sharp velocity profile. Through these 
standard tests, we confirmed that our code is capable to reproduce analytical results accurately and 
that all source terms other than the neutrino-matter interaction terms are correctly implemented. 

We then performed a practical core collapse simulation with which we could also confirm 
whether the above mentioned neutrino number conservation, especially in the trapping region, is 
satisfied adequately. We followed gravitational collapse, bounce and initial post bounce phase up to 
T p b ~ 50 ms of a 15 Mq star to make a detailed comparison with previous studies with Boltzmann 
neutrino transport. Regarding neutrino opacities, we currently employed baseline set where inelastic 
neutrino-electron scattering, thermal neutrino production via pair annihilation and nucleon-nucleon 
bremsstrahlung were included. We started the code comparison before core bounce. Important 
features such as the evolution of the central Y) and entropy as well as the neutrino trapping density 
showed nice agreement with the reference models. Next we made the code comparison after bounce 
till the stall of the bounce shock. At this transient phase, we checked that the Ml code can capture 
the overall evolution in the Boltzmann codes regarding the neutrino propagation from the opaque to 
the transparent region. Considering our insufficient resolution and the difference coordinates used, 
the neutrino luminosity and the energy spectrum of our code showed good agreements with the 
reference results. The peak luminosity {L Vetpea k = 3.9 x 10 53 erg s -1 ) showed almost the same value 
with the reference models. After the neutronization, the neutrino luminosities showed consistent 
values though there is a systematic lower shift in the electron-type neutrino luminosity. The RMS 
energies of the emergent neutrinos showed a similar level of the match with the reference results, 
which demonstrates the validity of our code. 

In the end, we shortly discuss for the future run (using Exa-scale platforms) in order to check 
a numerical convergence of our 3D results and to get a more closer match with the Boltzmann 
results in the much longer postbounce phase. As already mentioned, we were forced to adopt 
a low numerical resolution in this study (effective angular resolution behind the stalled shock 
~ 120 km is ~ 8km/120km ~3.8°) due to our limited computational resource. The code has been 
already tuned to get a high parallel efficiency (> 90 %, e.g., going from 2048 to 4096 cores) with 
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a very high performance efficiency (~ 32 %) which measures the ratio of the real performance of 
our code to the theoretical peak performance measured by the platform, Cray XC30 in our case. 
However, it still takes ~ 1.25 CPU days to follow ~ 1 ms postbounce in model 3DG15 by the peta- 
flops machine occupying ~ 10% of its total resource (20 48 processors'). Just only comparing the 
angular resolution < 2° of recent numerical studies (e.g., lHanke et al.1120131 : iMiiller fe Jankall2014l ; 


Takiwaki et ah 201 f). the employed resolution in this study is approximately two times coarser. 


Currently, the wall clock time per each timestep in the postbounce phase (the typical numerical 
time step At = 10 -7 s) is ~ 5 s if we use 4096 processors (i.e., we can follow the postbounce 
time of ^1.7ms per day, which we write as 1.7ms/day in the following). With this performance, 
we estimate how much computational resource and how many computational time are required for 
the twice high resolution model to follow ~ 200 ms after R 
expec ted to occur for low-mass progenitor stars in 3D (e.g., 


- 1 __ 

Takiwaki et al. 

(2014 

); 

Melson et al. 


(}201 5l)). Since our current resolution at the origin Ax ~ 480 m is marginally acceptable, we may 


just need to double the number of the numerical meshes in each nested box. 

If we are able to double the resolution with the fixed innermost mesh size and use 4096 x 2 3 ~ 
32,000 processors, we can also follow the same postbounce evolution ~1.7ms/day, or ~3.4ms/day 
occupying 4096x2 3 x 2 ~ 64,000 processors. Even if we can luckily take the latter case (~3.4ms/day), 
we still need approximately two months to follow ~ 200 ms after bounce occupying the 
processors. For more massive pi 


64,000 


1 -- 1 ___ A 

Marek et al. 

(2006); 

Muller et al. 

(2012b 

l); 

Nakamura et al. 


(2014bj)). This is apparently beyond the maximum computational time allocated to us in the K 
computer and surely needs Exa-scale platforms (in the next decade to come). 

Before the advent of these next-generation supercomputers, it is noted that we actually have 
lots of tasks to improve our code. We expect that we can still enhance the numerical efficiency 
especially when we get a convergent solution during the Newton-Raphson iteration in the im¬ 
plicit update. At present, the convergence becomes worse where the energy fluxes of neutrinos 
F^/H 11 become non-negligible (i.e., just above the neutrino sphere). We have confirmed that the 
neutrino-election scattering is one of the dominant factors which delays the convergence. With an 
eye towards actual application of this code in CCSN simulations, we plan to continue our code 
refinement, in which we not only need to employ a more el a borat e set of neutrino opacities (e.g .. 


Horowit j ( 2t)()dh 


Fischer et al 


Burrows et al 


J2013l l: 


( 2006 ;); Sumivoshi fc Ropke ( 2008 ): Martfnez-Pinedo et al. 


Bartl et al 


(20H 


( 2014) 47 with including the higher-order angular dependence of 


the reaction angle that was omitted for simplicity (e.g., Sec. I5TD . but also to find a more effi¬ 
cient algorithm to deal with the resulting (more complicated) Jacobi matrix in the implicit update. 
This study is only our very first step towards a more realistic coding, which is indispensable for 
quantitative study of stellar core-collapse and the explosion mechanism. 
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A. Neutrino Matter Interaction Terms 

In this appendix, we summarize neutrino matter interaction processes included in this study 
which are: absorption and emission process 

u e n -H- e~p, 
u e p <->• e + n, 
u e A -H- e~A', 

isoenergy scattering of neutrinos off nucleons and heavy nuclei 

un •H- un, (A4) 

up -H- up , (A5) 

uA -H- uA, (A6) 

inelastic neutrino electron scattering 

ue ue , (A7) 

thermal neutrino pair production and annihilation 

e _ e + uu, (A8) 

and nucleon-nucleon bremsstrahlung 

NN o NNuu. (A9) 


(Al) 

(A2) 

(A3) 


Four vector neutrino matter interaction term is given by summation of all these interaction 
terms as 


gfi _ ^na_|_ giso,fi _|_ gnes,fi _|_ gtp,/d _|_ gbrem,/x 


(A10) 


where 5 nae,/i , S lso ’^ ; 5 nes,A1 , 5 tp,/i and 5' brem,p are the four vector source terms of neutrino absorption 
and emission (nae), isoenergy scattering of neutrinos off nucleons and heavy nuclei (iso), inelastic 
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neutrino electron scattering (nes), thermal neutrino pair production and annihilation (tp) and 
nucleon-nucleon Bremsstrahlung (brern), respectively. 


We briefly summarize each interacting kernels, opacities and source terms in followi ng sub ¬ 
sections. For more detailed e xpressions and explanat ions, the reader is referred to iBruennl (119851 ) ; 
Hannestad fc Raffelt ( 1998 ): Rampp fc Janka (2QQ2J) . Derivation of the four vector source term 


from these quantities is given bv lshibata et ah ( 2011 1. 


A.l. Neutrino Absorption and Emission 


We consider that neutrino absorption and emission: by free neutrons v e n « e p, by free 
protons v e p <->• e + n and by heavy nuclei v e A <->• e~A'. The four vector source term S naB,fl is given 
by 


^nae,/i ^.nae 




v,e) 


47 T £ 3 

exp{/3(e — p u )} + 1 


0,D 


(All) 


where n nae is the opacity, f3 = l/ksT with &b the Boltzmann’s constant and ji Ve = —/rp e = 
/r e — p,p + fi n is the chemical potential of neutrinos in thermal equilibrium with matter. Here, pt e , 
[i n and p p are the chemical potentials of electrons, neutrons and protons, with including rest mass 
energy of each particle, respectively. For each reaction (z/ e n <->• e~p, v e p -H- e + n and v e A e~ A'), 
we first evaluate absorptivity 1/A [s -1 ] and then emissivity j [s^ 1 ]. 


Absorptivity for reaction u e n o e pis expressed as ()Bruennlll985l : iRampp fc Jankal 1200211 


1 

AT 


(he) 


(dv 4 " ^ C/A)Vnp 


7r 


l — F e -(e + Q) (e + Q)y/(e + Q ) 2 - 


rule 4 . 


(A12) 


where we employed gy = 1 and gA = 1-23 as form factors resulting from the virtual strong 
interaction processes. F x (e) = [1 + exp (f3(e — /ix))]” 1 is the Fermi distribution function of fermion 
x with energy e. Gp(= 8.957 x 10~ 44 MeV cm 3 ) is the Fermi constant. 

For reaction v e p -H- e + n, 


1 


f WF n ( 9 v + 3 g 2 A )ripn 

1 - F e + (e - Q) 

1 o, 



rriic 4 . 


for e > m e c 2 + Q^ A13 ^ 
otherwise, 


where Q = m n c 2 — m p c 2 is the rest mass energy difference of the neutron and proton. 
Finally for reaction v e A -H- e~A', 


1 

AT 


I W^A^N r N h nA 

1 - Eg- (e + Q') 

1 o, 



otherwise, 


(A 14 ) 
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where Q' = /i n — n P + A is the mass d ifference between the in itial and the final states through the 
reaction. By following iBruennl (|1985l l ; iRampp Sz Jankal (120021 ) , we employed A = 3 MeV and 


0 


N p = 


N h = 


for Z < 20 
Z - 20, for 20 < Z < 28 
8, for Z > 28, 

for N < 34 


6 


40 — N, for 34 < N < 40 


(A15) 


(A16) 


0, 


for N > 40, 


for the number of protons N„ and holes Nh in the dominant GT resonance for the electron capture. 


_ As for the values r] pn and r] np , we adopted ones proposed in IBruennl (|1985T ): IRampp fc Janka 

d2002l ) 


T]pn — 


Tlr 


exp 


l^p 


- 1 


^Inp — 


while we set 


exp (-/3(> n -n p - 


?Ipn — ^Ipi 
TJnp = 


- 1 


(A17) 

(A18) 


(A19) 

(A20) 


in the non-degeneracy regime where fi n — fi p — Q < 0.01 MeV is met. 


Once we evaluate the absorptivity, the emissivity j and the opacity k are obtained as fjBruenn 


1983 ) 


j e = e /3(e-/i„)_L 

A F 


(A21) 


and 


nae _ • i 
\e) “^+ Ae - 


(A22) 


A.2. Isoenergy Scattering of Neutrinos 

The four vector so urce term for isoe nergetic scattering of neutrinos off free nucleons and heavy 
nuclei is written as (Sh ibata et ah 2011) 


OlSO,/Z _ . iSO ttP 

b {e) - X{e) H (ey 


(A23) 
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To evaluate yjAj from the isoenergetic scattering kernel R lso (e, u>), we expand it into a Legendre 
series in term of u (here oj is the cosine of the scattering angle) up to the first order as 


v ’ ’ 2 2 ( £ ) ’ 


(A24) 


where and 4> 1 ^°’ 1 are the zeroth and first order of scattering kernels, respectively. After 

angular in tegratio n of the angular dependent source term with respect to ui, is expressed as 
below (Shibata et ah 2011) 


X(e) 


is° _ g 2 I ^,iso,0 _ (Jjiso,! 


' (0 


' 00 


(A25) 


For scattering process on free nucleon ( n/p ), zeroth and first order kernels become 




2t r Gi 


( £ ) (he) 3 2irh 


N \ 2 


N \ 2 


vnn uwy+3(h%y , 


w = UhVf -(«)-). 


( £ ) (he) 3 2irh 


(A26) 

(A27) 


Obviously, has a dimension MeV 2 s 1 and thus has a dimension s 1 . In above, N 

takes n (neutron) or p (proton) and hy and h% are defined as (Bruenn] 119851 ) 


h n — -- 
n V — 


hy = - - 2sin 2 (V, 


1 


h n A = -2«M, 

h A = \<]A- 


(A28) 

(A29) 

(A30) 

(A31) 


where 9 w is the Weinberg angle and we adopt the value sin 2 0w = 0.2325. By followin dMezzacappa Sz Bruenn 
(1993c); Rampp & Janka (2002), we evaluate Vnn as 


VNN = 

6v 

n N , -, 

(A32) 




£,n = 

3 

2 PE*' 

(A33) 

II 

^ 2 (3An N ) 2 ' 3 • 

2m jyc* 

(A34) 
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Next, for coherent scattering on heavy nuclei, zeroth and first order kernels become (jBruenn 


1985 ) 


= 

00 

^iso, 1 = 

00 


2tt G\ 
(he) 3 Anh A 
2?r G 2 f 
( he) 3 47t/i A 


(ac°- 

(ac%- 


( h A ~ z ) cly ) 


^ 2y- l + e~ 2 ^ 

^ 2 — 3t/ + 2t/ 2 — (2 + y)e~ 2y ' s j 


(A35) 

(A36) 


with = (h P V j A + h™ /A )/2, C^ /A = (h p y/A - h™ /A ), y = 4be 2 and b = 4.8 x 10 6 A 2 / 3 . Here n A , 
A and Z denote number density of heavy nuclei, the atomic number and the charge, respectively. 


A.3. Neutrino Electron Scattering 


Inelastic scattering of neutrinos off electrons plays important role to help neutrinos to escape 
more freely from centre as a consequence of down-scattering and it thus enhances deleptonization 
of central core . The source term , S^j is expressed in term of the collision integral B^ s 0 ^ as 
( Shibata et ah 2011 ) 


(e,n) 


= e3 / d^B^ Q) (u p + n, (A37) 

where U l is a unit normal four vector orthogonal to u p . The collision integral B^ s ^ along the 
propagation direction is expressed as 

B ?:h) = f £ ' 2 d£'dn' [f(e', fi') {1 - f(e, II)} iT es ’ in (e, ef,u) 

- {1 - f(e',Q')} f(e,n)R nes ’° u \£,e',u)] , (A38) 

where i? nes ’ in / out (e, e' ,u) is the angular dependent inward/outward scattering kernel and uj is the 
cosine of scattering angle, i.e., angle between 14 and Q'. With expanding the angular dependent 
inward/outward scattering kernel i? nes,m,/out (e, e', w) into a Legendre series with respect to co and 
taking up to the first order as 


D nes,in/out/ J \ 1 *nes,in/out,0 3 ,nes,in/out,1 

K (£,£ ,u) ~ 2 $ ( e,e0 + 2 

and with decomposing the neutrino distribution function after scattering into isotropic and non¬ 
isotropic parts as 



/(e',n , )«/ 0 (e , ) + / 1 ’ # *(e / K» 

final expression of the scattering integral is described as below 

= /( £ ,^) (A( £ ) ines + B (e),nes l v) + ^ej.nes + C {e),njv 


(A40) 


(A41) 
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In the above, we used the same notations used in Bruenn (1985) (see their Eqs.(A34)-(A39)) for 

°e),nes and C \e ),nes haVe three S P atial 


A9, , 

(£j,nes’ (e),nes 


C ( 0 £ ),ne S aIld C (e ),nes with tlle exception of Sj 


components (the zeroth component is deter mined fro m orthogonality = 0) and of -B^ nes 

a factor of 3 larger than that of Eq.(A35) in Bruennl (19850. They are explicitly written as 


is 


4° 

^ i (e),nes 


B 


_ 

(e),nes 


2vr 

(he) 3 

2ir 


e' 2 de' 

fi 


/ 0 ( £ ,) $ nes,in,° + (1 _ f( £ '))^ 
nes,in, 1 ^x-nes,out, 1' 


/>0 / /\\^nes,out,0 


r o 

°(e),nes 

C}’? 

(e),nes 


(he) 3 1 £> d£ 'f ^ W2) - 

2?r ^ AV Julies,in,0 


(W J f 


2vr 


where we employ 


/°(e) 


J (g) 
47T£ 3 ’ 


3/B\ 
_(fO 

47T£ 3 


(A42) 

(A43) 

(A44) 

(A45) 

(A46) 

(A47) 

(A48) 


for the isotropic and non-isotropic parts of the distribution functio n (see alsolShibata et al. ( 2011m. 
For an explicit evaluation for ( j ) “ s j™/ out ’ 0 / 1 ; we re f er the reader to Yueh &; Buchler (1976); Bruenn 
dl98, r Jl . 


Since coefficients (lA42h - (jA45[i do not depend on propagation direction 12, the final integral 
Ea. ()A37l) with imposing Eg. (I A4 111 becomes 


S 


nes ,/z 

(e) 


(■V" + Hfo) 4) jn „ 

+ (BaW* + L A.)) B l) 

+4 m 3 u“C° {e) ^ 


g (ernes’ 


(A49) 


here h )W = g^ u + u^u v is the projection operator. Since, <f> nes has a unit [cm 3 s 1 ], AjE nes) i?^ nes , 
C? \ and C}’? thus have a unit [s _1 l which is required to let the dimension of *S? e , s,Al to be 

(e),nes (e),nes l j i (e) 

[MeV 3 s -1 ]. 
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A.4. Thermal Pair Production and Annihilation of Neutrinos 


As for the thermal pair production/annihilation process of neutrinos, we take the same ap¬ 
proach as neutrino electron scattering process. The collision integral is expressed as 

B (?h) = f e' 2 de'd^[{l-j(e\^)}{l-f(e^)}R t ^ 0 (e,e'^) 

, (A50) 

where / denotes anti-neutrino distribution function and i? tp,pr0//ann (e, e', w) is the angular depen¬ 
dent production/annihilation kernel. We again expand the kernel i? tp,pr0//ann (e, e', w) up to the first 
order in uj as 

i . , „ , 

(A51) 


fltp.pro/ann, / ) „ I^p.pro/ann.O 3 tp.pro/ann,! 

\ i / 2 (£,£') 2 \f..e ) 


■ (e,e') 

and hnal expression of the scattering integral is thus described as below 

flgn, = /(*.«) (4),.p + B w,*pV) + c (p),.P + c ( ‘V- < A52 > 

which has exactly the same expression as that of neutrino electron scattering (Ea lA4ll) . Again 
coefficients have the same notations as in Bruenn (1 19851 ) (see their Eqs.(A43)-(A47)) and described 
as 

2vr 


4° 

A (<0,tp 


(he) 3 
2vr 




R 0 ’^ = _ 

( £ )4p (he) 3 


^(s),tp 

/WM _ 

U (e),tp 


Je’ 2 de'p*(e') , 


2t r 


2tt 


(^3 J e' 2 de'f^(e')^ 

For an explicit evaluation for <|)T’ pi °/ ann ’ 0 / 1 ) we re fer the reader to Bru enn (1985). The final 


(A53) 

(A54) 

(A55) 

(A56) 


expression of the source term is simply obtained by replacing coefficients in Eq. (IA49l) with 
those of thermal process, i.e., with Eas. (lA53l) - ()A56l) . 


A.5. Nucleon-nucleon Bremsstrahlung 


The collision integral for Nucleon-nucleon Bremsstrahlung has the same expression as 

that of thermal pair production/annihilation of neutrinos and, therefore, the same notations used 
in Sec. IA.4l can be directly applicable. is written as 


Bfc n) =J ^^[{1- f(J,Q')}{l - f(e,n)}R*^(e,6',u) 

-f(e',n')f(s,m hr ’ ann (£,£'^) 


(A57) 
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here ii br ’ pro,/ann (e, e', a;) is the production/annihilation kernel. As for the Bremsstrahlung process, 
we again expand the kernel i? br,pro/,ann (e, e', w) into a Legendre series, but take only the zeroth 
order term in w for simplicity. Then the kernel becomes 


^br.pro/ann^/^ 


1 t br,pro/ann,0 

2 ( £ ’ e ') 


(A58) 


Coefficients used in the final expression for the collision integral and the source term are simply 
evaluated by replacing <J? tp in Eqs. (|A53ft - (lA56l) with <f> br . 

Isotropic production kernel <J> b ^ ,pro ’ 0 can be evaluated by following manner. 


ibr,pro,0 


= e P ( £+£ ) ^2 m i n 

iGnn,pp,np 




(A59) 


In the above equation, indexes D and ND denote degenerate and non-degenerate limit of free 
nucleons, respectively, and (nn, pp, np) represent bremsstrahlung due to neutron-neutron, proton- 
proton and neutron-proton pair, respectively. e') and 0j^ D,, (e, s') are expressed as 


^b/Ob 6 ') = ( hc ) { 


\G 2 ( 2vr/ \ 4 


4 \ m 


mj 


x 


1 (47 T 2 +/5 2 (s + e') 2 ) 


OLi 


127T 9 /3 2 (e + e') 1 — e~P( £+£, '> Si 


2 he ( 3ir 2 Xjnb ) 
(he) 6 


,ND,it / \ n6 G 2 /2vr/\ 4 

K •(!,*) = &) T ( —) 


2m 


3/2 


n ll/2 pl/2( £ + £ ,) 


——— 2 x/TT^F+e 7 )' -G- X 2 n 2 . 


(A60) 
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In the above equations, (j)^J ND ’ 1 (£,£') have a dimension in [cm 3 s b ], / = 1 and G 2 = c(hc) 2 G 2 F = 
1.55 x 1CP 33 [MeV~ 2 cm 3 s -1 ]. We employed values for ctj, A,;, mj and S t as 



( 3 9a \ 

( x n \ 


( m n \ 


(G 

( nn ^ 

OLi = 

*9a 

■ X t =\ x p 

, mt = 

m p 

, Si = 

4 , for i = 

pp 


{ Wa / 

\ mm(X n ,Xp) y 


\ y/rn n m p ^ 


(G 

y np J 
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After deriving the production kernel, we can obtain the annihilation one by using a following 
relation 


i? br ’ ann ( £ , £ ',o;) = e^ £+£, )i? br ’ pro (e,e',w). 


(A63) 
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A.6. Summary of the Source Term and Mean Free Paths 

By combining all the source terms mentioned above, the final expression becomes 



+ V { (/O 1 ^ + H (e)) - 4 (hx + 4«VC ( “ )iX 


X Gnes,tp,br 


+ {H a (e) ufi + -^a M (£)) ^(e),X 


47T£ 3 


7 a 

' La °(e),X 


where is written by 


■ft 


47 T£ 3 

exp{/3(e - //„)} + 1' 


(A64) 


(A65) 


In the end, we plot inverse mean free paths of each process for reference. Assumed hydrody- 
namical profiles are p = 10 10 g cm^ 3 , T=0. 638 MeV and Y e = 0.43 (Fig. H71) . which corresponds to 
the central profile of “s!5s7b2” in Woosle y fc Weaver (1995), p = 1 0 11 g cm~ 3 , T=1.382 MeV and 


Y e = 0.4 (Fig. fTH]) . which is the same condition used in Fig. 36 in iBruerm] (1985) and p = 10 13 g 


cm 


-3 


T=3 MeV and Y P = 0. 25 (Fi g. fT9l) . which is a typical profile after neutrino trapping during 
collapse phase. As in Bruennl (1985), we assumed the final state occupancy of the neutrinos is zero. 
More explicitly, we plot 1/A m f p [cm -1 ] expressed by 


A 


-l 

mfp 



V 1S ° 

X (ep 

A9, , 

(ernes’ 

4° 

We),tp’ 

4° 

(e),br’ 


(neutrino absorption and emission) 
(isoenergy scattering) 

(neutrino electron scattering) 
(thermal pair production and annihilation) 
(nucleon — nucleon bremsstrahlung). 


(A66) 


Since isoenergy scattering rate has the same value for all neutrino species, we only plotted for 
electron type neutrino. Note that, the inverse mean free path for bremsstrahlung is multiplied by 
10 15 (Fig. fT71) . 10 10 (Fig. fT8l) and 10 5 (Fig. fT9l) and, for thermal process, it is multiplied by 10 10 

(Fig. mm- 
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1, i/A^i/A 5, ^ e n^e p(^ e p^e + n) 

2, vn^vn 6, ve^ve 

3, i/p^i/p 7, i/i/NN^NN (xIO 15 ) 

4, i/A?^e - A' 8, i/i/^e~e + (xIO 10 ) 



Fig. 17.— Energy dependence of inverse mean free paths for each process labeled at top part. 
Employed hydro dynamical background is p = 10 10 g cm -3 , T=0.638 MeV and Y e = 0.43 
and we assumed t he fi nal s tate o ccupancy of the neutrinos is zero. With using the EOS of 


Latti mer & Douglas Swesty (]199l|) with compressibility parameter K = 220 MeV, several rep¬ 
resentative thermodynamical quantities become s = 0.71 k-Q baryon -1 , A = 65.5 and Z = 28.2. 
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1, i/A^i/A 5, i^n^e p(i' e p? ± e + n) 

2 , i/r\^i/n 6 , ve^ve 

3, i/p^z/p 7, j/j/NNMNN (xIO 10 ) 

4, vk^e~K 8, w^e~e + (xIO 10 ) 



Fig. 18.— Same as FigfTTl but with p = 10 11 g cm" 3 , T=1.382 MeV and Y e = 0.4. This hydro- 
dynamical profile returns s = 1.19 k b baryon” 1 , A = 73.9 and Z = 30.5. Since N = A— Z > 40, 
neutrino absorption on heavy nuclei is blocked with our current electron capture rate fEa lA16l) . 
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1, i/A^i/A 5, i^ e n? i e - p(i7 e p? i e + n) 

2 , i/r\^i/n 6 , ve^ve 

3, l/p^up 7, j/i/NNMNN (xIO 5 ) 

4, i/A?±e - A' 8, w^e~e + (xIO 10 ) 



Fig. 19.— Same as Fig(17]but with p = 10 13 g cm” 3 , T=3 MeV and Y e = 0.25. This hydrody- 
namical profile returns s = 1.25 k-Q baryon” 1 , A = 167.7 and Z = 49.3. Since N = A — Z > 40, 
neutrino absorption on heavy nuclei is blocked with our current electron capture rate lEa lA16l) . 
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B. Comparison in ID Spherical Coordinate 


In this appendix, we perform a core collapse simulation using ID spherical symmetric Ml 
neutrino transport code. The aim of this section is to test the local neutrino-matter interaction 
terms implemented in our main 3D Cartesian grid based radiation-hydrodynamics code. Since 
our 3D multi-energy neutrino radiation-hydrodynamics code with higher resolution than the one 
used in Sec. [5] is still computationally demanding, it is not easy to check the spatial resolution 
dependence on the practical core collapse simulation. Therefore, we can not assess whether the 
differences between our (pseudo-) ID results and previous ID ones come from the spatial advection 
terms or from the neutrino-matter interaction terms. To expel the ambiguities coming from the 
spatial advection term, we develop a new spherically symmetric Ml code in the Newtonian limit. 
The basis of the new ID code is the same as our main 3D Cartesian grid based code. We solv e 
Eg. (1321) with the same neutrino opacities used in model “sl5Nso_ld.b” in i Buras et al. (2006b|). 
The spatial advection t e rm S a ^, s in ID sphe rical symmetry is evaluated in a similar way as in 


Nakamura et al.l ( 2014b| ); Horiuchi et al- (2 0141 ) and the advection terms in energy space S a d v ,e are 
also slightly modified for spheri cal p olar c oordinate. To take the same approach as the recent Ml 
neutrino transport codes (O’Connor 2014; Just et al. 20151), we solve the advection terms in energy 
space S a dv, e explicitly in time. 

In Fig l20l we plot our results (ID—Sph) for the neutrino luminosity, mean energy, shock 
radius and the W profi le together with reference values of VERTEX-PROMETHEUS (model 


; ‘s!5Nso_ld.b” in Burasetal 


2006b| ) (VX—N15) except the Y e profile. In the Y e profile, we also 


plot results obtained from a model with solving S ac j v e implicitly in time, i.e., as the same as our 
main 3D code. Our results show a good agreement with VX—N15. The maximum deviations are 
~ 7x 10 51 erg s -1 in the luminosity and ~ 0.9 MeV in the mean energy. The mean energies, (£y e ) 
and (sp e ), show slightly higher values than the reference ones, they are, however, within the error 
bounds reported in previous code comparisons (ILiebendorfer et al.l 120051 : iMiiller et al.ll201(l ). The 
neutrino luminosities show quite consistent behavior especially in late phase T p b > 150 ms, when 
the neutrino heating becomes more important. Contrary to our main results reported in Sec. 15.21 
the Y e trough observed at R r\_; 10 6 cm (see, Fig. fT3l) disappears in the explicit model (solid lines), 
while the Y e trough still exists in the implicit model (dotted). These two models support the idea 
to evaluate both S a d v , s and S a d v , e at the same time slice. We are going to examine this point in 
our main 3D Cartesian based code in the near future. From this test, we confirm that the neutrino 
matter interaction terms are correctly implemented. 
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T pb (ms) 



Fig. 20. Top left panel: Time evolution of the neutrino luminosities measured at 400km. Results 
for our ID-spherical Ml transport scheme (ID—Sph) and for VERTEX-PROMETHEUS (VX—N15) 
are denoted by red and green lines, respectively. Top right: Same as the top left one but for the 
mean neutrino energies. Bottom left: Time evolution of the shock radius. Bottom right: Radial 
profile of the electron fraction Y e at different time slices. Solid and dotted lines are for explicit and 
implicit models, respectively. 
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